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Abstract 

The  direct  and  inverse  scattering  problem  is  solved  for  the  elliptic  Sinh- 
Gordon  equation.  It  is  shown  that  the  inverse  scattering  transform  may 
be  useful  in  the  analysis  of  localized  singular  solutions.  As  an  example,  a 
cylindrically  symmetric  singular  solution  is  discussed  in  detail. 
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1  Introduction 

The  elliptic  Sinh-Gordon  equation 

ur*  +  uvv  =  A2  sinh  u  (1) 

appears  in  plasma  physics  [1-5]  as  a  two  dimensional  model  equation  describ¬ 
ing  a  sj'stem  of  interacting  charged  particles.  Depending  on  the  sign  of  A2 
we  can  consider  both  positive  and  negative  temperature  states  of  thermal 
equilibrium  [4],  corresponding  to  essentially  different  solutions  of  eq.  (1). 

Negative  temperature  states  (A2  <  0)  have  been  studied  extensively  [3-5], 
and  both  numerical  and  analytical  solutions  have  been  reported,  exhibiting 
stable  nonuniform  distribution  of  the  potential  and  charge  density  within  a 
bounded  region. 

On  the  other  hand,  little  attention  has  been  paid  so  far  to  the  positive 
temperature  states  (A2  >  0),  perhaps  because  of  the  fact  that  the  only  pos¬ 
sible  regular  solution  (if  it  exists)  is  trivial,  i.e.  corresponds  to  the  uniform 
distribution  of  charge  density. 

The  case  A2  >  0  becomes  nontrivial  if  we  place  a  fixed  charge  (or  a 
number  of  charges)  into  the  medium  described  by  eq.  (1),  i.e.  if  we  allow 
the  potential  to  be  singular  at  some  point(s)  within  the  region  of  interest. 

In  this  paper  we  confine  our  attention  to  this  latter  case,  and  without  loss 
of  generality  we  put  A2  =  1,  by  rescaling  eq.  (1)  to  the  dimensionless  form 

u**  +  Uyy  =  sinhu  .  (2) 

A  linearized  version  of  eq.  (2) 

u*x  -f  Uyy  -  u  =  0  (3) 
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can  be  easily  solved  using  the  Green  function  formalism  [6].  Indeed,  for  a 
unit  point  source  at  f  and  the  boundary  conditions  u  — »  0  as  |f]  — ►  oo,  the 
Green  function  is  given  by 


G(r,fO  =  /M|r-f'|)  ,  (4) 

where  A’o(r)  denotes  the  modified  Bessel  function  of  the  second  kind  and 
order  0,  r,r  '  are  vectors  with  the  components  (x,y),(x',  y',  ),  respectively. 

Placing  a  point  charge  of  strength  A  at  the  origin  we  find  a  cylindrically 
symmetric  solution 

u(r)  =  AI<0(r ),  r  =|  r  |=  y/x*  +  y2  ,  (5) 

which  satisfies  eq.  (3)  everywhere,  except  for  the  singularity. 

On  the  other  hand,  in  the  strong  nonlinearity  limit  we  can  expect  u  -+  oo 
as  r  — ►  0,  and  the  solution  of  (2)  tending  to  that  of  the  Liouville  equation 


uXT  4  uvv  = 


(6) 


The  general  solution  of  eq.  (6)  is  given  in  terms  of  two  arbitrary  functions 
[7].  Imposing  cylindrical  symmetry  we  can  find  a  class  of  solutions  singular 
at  the  origin,  with  the  leading  terms  given  by 


u  =  -alnr  4  /?  4  0{r)  ,  (7) 

where  a  e  (0,2). 

Unfortunately,  in  contrast  to  eqs.  (3)  and  (6),  a  closed-form  solution  of 
eq.  (2)  having  cylindrical  symmetry  is  not  known.  Naturally,  one  can  use 
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numerical  methods  to  integrate  (2),  starting  e.g.  from  the  asymptotic  expres¬ 
sion  (5).  In  this  paper,  however,  we  apply  the  inverse  scattering  transform 
(1ST)  [8-11]  in  order  to  reconstruct  the  solution  of  (2)  for  arbitrary  f. 

It  should  be  noted  here,  that  the  elliptic  problems  are,  in  general,  well- 
posed  when  the  boundary  conditions  are  imposed  along  a  closed  curve  sur¬ 
rounding  the  region  of  interest  [6].  On  the  other  hand,  the  initial  value 
problem,  which  is  a  natural  choice  for  the  hyperbolic  equation,  may  be  un¬ 
stable  (ill-conditioned)  when  applied  to  the  elliptic  case.  In  this  context,  wfe 
could  expect  possible  difficulties  in  the  application  of  the  1ST  to  an  elliptic 
problem.  In  fact,  the  1ST  has  always  before  been  used  only  for  solving  evo¬ 
lution  equations  of  the  hyperbolic  type,  and  to  the  authors  knowledge  it  has 
not  been  adapted  to  the  elliptic  case. 

Therefore,  the  main  aim  of  this  paper  is  to  study  the  applicability  and 
usefulness  of  the  1ST  formalism  to  the  analysis  of  elliptic  problems.  In  this 
connection,  in  Section  2  we  consider  the  direct  scattering  problem  in  a  rather 
general  case,  without  specifying  details  of  the  solution.  In  particular,  the  an¬ 
alytic  properties  of  the  scattering  data  are  determined  by  following  closely 
the  methods  of  Refs.  [9,11].  The  y-dependence  of  the  scattering  data  is 
discussed  in  Section  3,  while  Section  4  deals  with  the  inverse  scattering  prob¬ 
lem.  In  Section  5  we  consider  a  cylindrically  symmetric  solution  of  eq.  (2), 
having  singularity  at  the  origin  and  vanishing  exponentially  as  r  — +  oo.  Sec¬ 
tion  6  contains  concluding  remarks,  and  we  shall  also  point  out  some  open 
problems. 


2  Analytical  Properties  of  the  Scattering  Data 

Following  [8,10],  let  us  consider  the  linear  eigenvalue  problem 


(C  coshu\  (  sinhtA 


(8a) 


(sinhu\  .(C  coshu\ 

?+-— r+,U-— r”  (86) 

where  p+  =  \(iux  -f  uv)  and  £  is  the  spectral  parameter.  Take  the  y- 
dependence  of  the  eigenfunctions  t>i,  t;2  to  be 


vi,v 


(C  coshtA 

2  +  “8C~j 


v\  —  i 


sinh 

“sT; 


V2  , 


(9a) 


(sinhtA  (C  coshu\ 

P+  +  — Jv,-^+— .  (96) 

Then  it  can  be  verified  by  cross-differentiation  that  (8a, b)  and  (9a, b)  are 
compatible  if:  (i)  u  satisfies  eq.  (2)  and  (ii)  £  is  independent  of  y. 

For  C  real  and  u  tending  to  zero  sufficiently  rapidly  as  x  — ♦  ±00  we  define 
the  solutions  of  (8a, b)  with  the  asymptotic  form: 
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where  k(()  =  §  - 

It  can  be  shown  that  if 


is  a  solution  of  (8a, b),  then 
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u2(i,  -C} 
-C) 


is  also  a  solution. 

In  particular,  we  have 


4>{xX) t'kx  as  x  — »  —  oo 


(12a) 


V>(x,0-*  o  e~tkz  35  1  +°°  • 

Since  xj>  and  \j>  are  linearly  independent,  we  can  write  for  £  real 


(126) 
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(13a) 


^(x,C)  =  6(C)t/>(^,C)  ~  5(C)t6(i,C)  > 


where 


«(08(  0  +  6(C)6(C)  =  1 

Consequently,  it  follows  from  (11)  that 


(136) 


8(C)*«(-0  , 


(15a) 


6(0  =  6(-C)  • 


(156) 


One  can  also  show  that 
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“(jt)  .  a«<0 

\4C/ 

=4-(D  ,  (164) 

\4C/ 

and  if  u(x,—  y)  =  u(x,y)  then 

a(y,  0  =  n'(-y,  C)  ,  (Ha) 

My,  0  =  -»•(-».  H  .  (176) 


the  leLSt  two  relations  being  particularly  useful  when  the  solution  of  eq.  (2) 
is  symmetric  with  respect  to  the  t/-axis. 

In  order  to  determine  the  analytical  properties  of  the  eigenfunctions  <f>,  rp 
we  replace  the  differential  equation  (8)  satisfying  the  boundary  conditions 
(10)  by  an  equivalent  integral  equation  [9,11] 

*1(I,C)e“<*<)=l  +  f  N>(x,z, {)#,(*,  ()«““■•<>&  ,  (ISa) 

J —  oo 

*(*,C)e“(t'°  =  /'  A',(x,2,C)*(z,C)e“''’<><l2  ,  (1  St) 

J  —  oo 

where 

JV.(*,«,C)  =  ?(x)  , 


1 


6 


9(x)  =  M*)-^|HW,  rW  =  _(Ml)  +  !i£^£) 


<*(*>()  =  *(0*  -  J_°_  sink1  '~-is  . 

Extending  d>i(x,£)  into  the  upper  half  plane  ((  =  £  +  ir),T)  >  0)  one  can 
show  that  for  |  £  |>  | 

|  £)e‘°'(r,(‘)  |<  e7^coshV+(x)  ,  (19°) 

and 

|  &(x,C)e,o(x-°  |<  e7(x)sinh  V+(x)  ,  (196) 

where 

H*)  =  '—■dz 

T/+(i)  =  J  oo[  I  P+(z)  t  +\  I  sinhu(z)  |J  dz  . 

For  I  C  l<  I  we  should  consider  an  alternative  form  of  the  eigenvalue 
problem  [9].  Proceeding  as  before,  we  find  for  i)  >  0 

|  <j> i(x,Qe'kx  |<  e7^  1 1  cosh  j  |  cosh  V_(x)+‘|  sinh  ^  J  sinhVL(x)|  , 

(20a) 

|  fa(x,()e,k*  |<  e7^  1 1  cosh  ^  |  sinh  Vl(x)+  |  sinh  ^  |  cosh  VL(x)j  ,  (206) 
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where  V_(x)  =  /f^  {|  P-(z)  |  |  sinhu(z)  ||  dz  and  p.  =  \  ( iux  -  uv). 

Differentiating  4>(x,()e,kx  with  respect  to  £  we  can  see  that  the  deriva¬ 
tive  exists  for  tj  >  0.  Thus,  we  can  summarize  the  analytical  properties  of 
eigenfunctions  4>,  t/>  and  scattering  coefficients  a,  6: 

Theorem  -  If  the  following  integrals 


r+oo 

r+oo 

r+oo 

/  1  1  dx, 

/  |  u„  |  dx, 

/ 

J  —  OO  j 

—  OO 

/  — OO 

are  bounded,  then  4>(x,^)etkx,il>(x,()e~>kx  and  a(()  are  analytic  in  the  upper 
half  plane  (tj  >  0).  For  t)  =  0,  the  above  functions,  as  well  as  6(()  are  at 
least  bounded. 

For  (  — ►  oo  in  the  upper  half  plane  one  can  find  that 


*(x,Oe“'  = 


(21a) 


and 


(216) 


Similarly,  for  ( 


<.(0  =  140(1 

0  in  the  upper  half  plane  we  have 


«*,C)e** 


cosh  j 
i  sinh  j 


+■0(0  . 


<H*,  O''*’ 


—i  sinh  j 
cosh  j 


+  0(0  , 


(21c) 


(22a) 

(226) 
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c(0  =  i  +  0(0  . 


(22  c) 


3  The  y-dependence  of  the  scattering  data 

Eqs.  (9a, b)  can  be  written  in  a  matrix  form  as 


where 


vv 


A  B 

C  -A  V  ’ 


(23) 


C(x,y, 0  =■  i  •  (24) 

Since  u(x,y)  — ♦  0  as  |  x  |— *■  oo,  we  have 


lim 

|r|—oo 


lim  B(x,y,()  =  lim  C(x,y,C)  =  0  .  (25) 

|r|— *oo  |r|— »oo 

Eq.(23)  is  satisfied  by  the  functions  <fi(z,y,()eAov  and  4>(x,y,()e~AoV,  thus 
the  corresponding  equations  for  <f>  and  <t>  can  be  written  as 

*•-[*0*  -A-A.]*  <26‘> 
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and 


<f>v  = 


A  ■+■  Ao  B 
C  — A  +  Aq 


<t> 


(2  66) 


Taking  the  limit  x  — *  -foo  in  (26a, b),  we  find  the  differential  equations  for 
the  scattering  coefficients: 


av(y>  0  =  0  , 

&v(y,C)  =  2A0b(y,()  , 
av(y,0  =  0  , 

My,  0  =  -2A0b{y,O  ■ 

Thus 

o(y,C)  =  a(0,C)  , 

i-(V.C)  =  i(0,C)e-2A»(<>»  , 

a(y,()  ~  a(0,C)  > 
i(y,O  =  5(0,c)e2^)t' . 


(27a) 
(276) 
(27c) 
(27  d) 

(28  a) 
(2S6) 
(28c) 
(28  d) 
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4  The  Inverse  Scattering  Problem 

Following  [9]  we  assume  rp(x ,£)  to  be  given  by 


^(*,0 

where 


0 

1 


JC(x,s)  +  —M(x)L(x,s) 


eik,ds  , 


(29) 


cosh  j  i  sinh  | 
tsinhj  —  cosh  j 

and  K,  L  are  independent  of  £. 

Substituting  (29)  into  the  eigenvalue  problem  (8a, b)  and  taking  the  Fourier 
transform  we  obtain  the  following  differential  equations  for  K.  and  L : 


1  ufi) 

[Idx  —  er3d,  +  tCT2p+(i)]AC(i,s) -f  -cr2  sinh -L(x,s)  =  0  ,  (31a) 


[Idx  -  <73dt  +  ia2p_(i)]L(z,s)  +  ^a2sinh  ^^AC(i,s)  =  0  ,  (316) 

subject  to  the  boundary  conditions: 

\nnK(x,s)  =  0  ,  (32a) 

limi(i,i)  =  0  ,  (326) 

=  \p+{*)  ,  (32c) 
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£i(x,z)  =  -isinh^r  ,  (32c?) 

were  ct:,  02,  03,  denote  the  standard  Pauli  matrices. 

We  can  see  that  the  solution  for  K.  and  L  exists  and  is  unique  (except  at 
the  singularities  of  u(x)). 

In  order  to  derive  the  inverse  scattering  equations  we  consider  the  follow¬ 
ing  integral  in  the  complex  £  plane  (for  £  below  C ): 


Jo 


dC<j>(x,C)tik{C}x 

(C  -  CMC) 


(33) 


where  C  denotes  the  contour  extending  from  —  00  +  i0+  to  0",  then  from  0+ 
to  -foo  -f  i0+,  and  passing  over  all  zeros  of  a(('). 

For  u  on  compact  support  the  Jost  functions  as  well  as  a(C),  HO  are 
analytic  everywhere  (except  for  £  =  0  and  £  =  00),  and  we  can  express 
<f>{x,£')  in  terms  of  r/>(x,  (')  and  ^(x,£') 


f  dC'*(z,ne*ia*  r 

Jc  (C'-C)a(C')  JcC-C*{'() 

dC 


)* 


+ jcT~rrP(0)^xxyHC)x , 


C'-C 

where  p(Q  =  &(C)/«(0* 

Using  (21)  and  applying  the  Cauchy  theorem  we  find 


(34) 


V»(x,C)e 


»*(<)*  _ 


+ 


1 


■  US. 1 


p( C')0(x,Oe*l°*  •  (35) 


2*,JcC-C- 

If  u  is  not  on  compact  support,  the  above  contour  integral  can  be  replaced 
by  an  integral  along  the  real  axis  plus  all  contributions  from  the  poles  of  /?(£). 

Note,  that  V>(x,£)  can  be  also  expressed  in  terms  of  K  and  L ,  by  using  (11) 
and  (29).  Thus,  substituting  (29)  into  (35)  and  taking  the  Fourier  transform 
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we  obtain  the  inverse  scattering  equations  of  the  Gelfand-Levitan-Marchenko 
(GLM)  type  for  z  >  x  : 


ic2fC(x,  z)  -f 


0 

1 


F(°)(x  +  z)+ 


J°°  ds  [£(x,  s)F(0\s  +  z)  +  £(x,  s)F(1)(s  +  *)]  =  0  ,  (36a) 


icr7C(x,z)  -f 


0 

1 


FV(x  +  z)+ 


J~  ds  [/C(x,  s)F(1)(s  +  z)  +  £(x ,  s)Fw(s  +  z)\  =  0  ,  (36 6) 


where 


C(xys)  =  M(x)L(x,s)  , 


(37) 


(38) 


Once  K.  and  £  are  found  from  (36a, b),  the  solution  u  and  its  y-derivative 
uv  can  be  recovered  using  (32c, d)  and  (37) 


tanh^  =  TiA(l’l) 


2  i  +  £2(x,x) 


(39a) 


uv  =  8£i(x,x)  —  iux  .  (396) 

5  An  Example 

According  to  the  above,  if  the  solution  u  vanishes  as  x  — ►  ±oo,  then  we 
can  reconstruct  it  from  the  scattering  data  at  any  value  of  y  where  u  and  its 
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derivatives  are  nonsingular  for  all  x.  But  how  do  we  determine  the  scattering 
data?  And  how  do  we  determine  its  evolution  across  any  value  of  y  for  which 
u  has  a  singularity?  Before  we  can  answer  these  key  questions,  let  us  first 
study  the  simplest  nontrivial  example  of  a  solution  u  which 

(i)  is  cylindrically  symmetric, 

(ii)  is  singular  at  the  origin, 

(iii)  vanishes  sufficiently  rapidly  as  r  — ►  oo. 

A  closed  form  of  the  above  solution  is  not  known;  however,  for  r  — »  oo 
we  can  easily  find  that  the  solution  approaches  the  linear  limit  (5): 

u  ->  AI<0{r)  ,  (40) 

where  I\  o(r)  denotes  the  modified  Bessel  function  and  A  is  a  real  constant. 

If  u,  ux  and  uv  are  infinitesimal,  then  the  scattering  coefficients  can  be  cal¬ 
culated  as  standard  (linear)  Fourier  transforms  [9,11]  of  u  and  its  derivatives. 
In  particular,  from  (18a,b)  we  find  for  u  ■<  1 

a(y,C)  “+  1  »  (41a) 

Ki/.O  -  /+“  [m*,v)  -  ■  («») 

Expressing  u  and  p+  in  Cartesian  coordinates  for  r  — ►  oo  we  have 

u(x,y)  =  AI<o  [yjx7  +  yA  ,  (42a) 
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X 


(426) 


«.(*,»)  =  ft  (v/?TP)  , 


«*(*,!/)  =  -A- 


rl< 


i  (\/*J  + !/’) 


(42c) 


v/z2  +  yi 

The  Fourier  transforms  of  the  above  functions  can  be  calculated  by  using  the 
method  of  Ref.  [12].  For  k  real  and  y  >  0  we  find 


[+°°  Ko  (y/x7  +  v2)  e'2ikxdx  =  -=L  ...  . 
J-<*>  J  +  (2  ky 


,-vvW 


r  K'  = vww , 

•'-«>  V*  +  jr  1/ 

where  yjl  -f  (2 A')2  is  assumed  to  be  positive. 


(43a) 


(436) 


Substituting  (42)  into  (41b)  and  using  (43a, b)  we  obtain  for  y  >  0,  (  real 
and  positive  (f  >  0) 


Hs/.C)=->lfe-''v,’+(’*>!  ,  (44a) 

while 


£(!/,<)  =  0.  (446) 

Since  the  scattering  data  are  not  analytic  for  y  =  0  or  for  (  =  0,  we  should 
consider  separately  the  case  y  >  0  and  y  <  0  as  well  as  £  >  0  and  £  <  0. 

For  y  >  0  and  (  on  the  negative  real  axis  ((  <  0)  we  have 

&(y,O  =  0  *  (45a) 
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Hv,C)  = ■  («6) 

Similarly  we  can  find  for  y  <  0  and  (  >  0: 

Hu,  0-0,  (46a) 

=  ,  (466) 

while  for  y  <  0  and  £  <  0 

Hv,<  )  «  x|e*VwI>?  ,  (47a) 

Ky,  0  =  0.  (476) 

On  the  other  hand,  it  follows  from  (28, c)  and  (41a)  that 

«(y»  0  =  «(y»  0  =  1  (48) 

for  any  y  ^  0. 

Note,  that  the  expressions  (44)-(48)  for  the  scattering  data  are  in  agree¬ 
ment  with  both  (28)  and  (15a, b), (17a, b).  Moreover,  since  we  can  make  u  as 
small  as  we  wish  by  letting  y  —*  ±oo,  we  can  conclude  that  the  above  results 
are  exact  and  valid  for  arbitrary  y  (except  for  y  =  0). 

Having  determined  the  scattering  data,  we  are  in  a  position  to  consider 
the  inverse  scattering  problem.  Since  a(y,()  =  1,  there  are  no  bound  states, 
and  contour  integral  in  (38)  can  be  replaced  by  an  integral  along  the  real 
axis: 
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F^=hC«{hYp«)e‘m‘  *  (49) 

where  p(0  =  6({)/a(0  and  k(£)  =  f 

Substituting  (44a),  (45a)  and  (48)  into  (49)  and  performing  integration 
we  find  for  y  >  0: 


'"«~7  (’♦*)-£ 


*■  V  f  +  »’ 


iA„  /  //e\2 


^-7*  V  I  +  *’  ' 


(50a) 


(505) 


vnv(*)+*’ 

^’w-s(»-2)— 77^r- r 


f  +y2 


(50c) 


For  y  —*  oo,  p(£)  becomes  infinitesimal  and  we  find  from  GLM  equations 
(36a, b): 


u  =  8t£i(x,x)  =  8iF^(2x)  =  AKo  [\Jx2  +  y2 )  , 


(51a) 


u„  =  8£i(x,x)  —  tit*  =  8F(0^(2x)  —  tu* 
-  Ay 

On  the  other  hand,  for  y  — ♦  0  we  note  that  the  limit 


lim  F(n)(z) 

v-.o+  v  ' 


(516) 


exists.  Moreover,  similar  calculations  performed  for  y  <  0  show  that 
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(52) 


lira  F<n\z)  =  lim  Fw0)  , 

y— *0“  y— ♦0+ 

in  spite  of  the  fact  that  the  scattering  data  are  discontinuous  for  y  =  0. 
Thus,  for  y  =  0,  we  have  simply 


F“»«  =  CKiC-)  , 

(53a) 

Fm{z)  =  CK 0(|)  , 

(536) 

F«\z)  =  CKi(2t)  , 

(53c) 

where  C  =  —iA/8,  and  it  remains  only  a  technical  problem  of  solving 
the  GLM  equations  (36a, b)  with  relatively  simple  expressions  (53a,b,c)  for 
F^{z). 

6  Conclusions 

In  this  paper  we  have  discussed  the  direct  and  inverse  problem  associated 
with  the  elliptic  Sinh-Gordon  equation.  In  particular,  attention  has  been  paid 
to  the  simplest  case  of  a  cylindrically  symmetric  singular  solution,  for  which 
we  have  been  able  to  derive  exact  analytical  expressions  for  the  scattering 
data  (44)-(48).  The  last  step  (i.e.  solving  the  GLM  equations)  cannot  be 
done  analytically,  however  due  to  simple  (and  exact)  expressions  for  the 
kernels  F^n\z)  the  problem  is  tractable  by  numerical  means  and  allows  us 
to  find  effectively  the  potential  u  and  its  derivatives  ux,uy. 
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Generally  speaking,  we  have  found  the  1ST  approach  to  be  surprisingly 
effective,  in  spite  of  a  rather  nonstandard  application  to  the  elliptic  problem. 
However,  some  questions  remain  open: 

(i)  Both  numerical  results  and  the  WKB  analysis  [13]  show  that  the  am¬ 
plitude  A  (see  eq.  (40))  is  bounded  by  Amax  =  4/7 r.  For  A  >  Amax  the 
solution  seems  to  enter  a  new  class  which  has  point  singularity  surrounded 
by  a  singular  ring. 

(ii)  Generalization  to  the  case  of  two  (or  many)  charges  is  nontrivial.  Pre¬ 
liminary  results  indicate  that  the  apparent  asymptotic  charge  distribution 
differs  significantly  from  that  of  the  single  charge. 

We  will  address  the  above  problems  in  a  future  publication. 
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Abstract 

An  analysis  of  the  initial  value  problem  for  the  planar  magnetron  reveals 
that  a  coherent  structure  (convective  cell)  can  be  created  from  arbitrary 
initial  conditions  and  that  these  structures  will  grown  linearly  in  time.  This 
is  proposed  as  the  explanation  of  the  origin  rf  the  convective  cell  formation 
seen  in  recent  numerical  Simula 


1  Introduction 

Although  there  has  been  a  very  long  history  of  the  study  of  the  eigenmodel 
structure  of  a  planar  magnetron1-5  as  well  as  also  cylindrical  magnetrons6-8, 
to  our  knowledge,  no  study  has  been  made  of  the  general  initial  value  problem 
for  infinitestimal  perturbations  for  these  devices.  We  do  find  that  the  general 
features  of  such  a  treatment  for  low  density  nonneutral  plasmas  has  been 
discussed9-11  and  is  well  known.  In  the  low  density  case,  it  is  found  that  any 
infinitestimal  disturbance  in  the  potential  will  tend  to  decay  algebraically  in 
time9-11. 

These  results  have  been  frequently  quoted  by  others12  as  sufficient  rea¬ 
son  for  ignoring  the  continuous  spectrum  also  in  the  high  density  region  of 
these  devices.  While  that  may  be  true  for  the  electrostatic  potential,  we  shall 
demonstrate  here  with  a  simple  asymptotic  expansion  that  such  is  not  so  for 
the  fluid  motion.  Indeed  we  find  that  a  very  important  generic  algebraic 
instability  always  exists  in  these  devices.  And  it  can  also  occur  at  low  densi¬ 
ties,  too,  although  reduced  in  magnitude  by  a  factor  of  u.,2/fl2.  Furthermore, 
this  result  will  demonstrate  that  "coherent  structures13”  can  arise  from  the 
plasma  fluid  equations.  The  generic  features  which  allow  them  to  exist  in 
this  case  are  the  shear  flow  and  the  wave-particle  resonance. 

The  simplest  example  of  a  solution  of  a  general  initial  value  problem 
for  infinitestimal  perturbations  of  a  fluid  flow  problem  has  been  given  by 
Case14'15.  As  to  details,  we  refer  the  reader  to  Refs.  [9-11,14,15].  Here  we 
shall  briefly  outline  the  solution  of  the  initial  value  problem  for  perturbation 
of  the  planar  magnetron  (cylindrical  case  will  be  essentially  the  same  except 
for  changes  due  to  the  cylindrical  coordinates  and  curvature  effects),  detailing 
only  the  important  differences  from  the  standard9-11,14,15  case. 
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2  The  Initial  Value  Problem  for 

Infintestimal  Perturbations  of  a  Planar 
Magnetron 

We  shall  model  the  planar  magnetron  with  the  cold-fluid  plasma  equations 
using  smooth  bore  boundary  conditions.  First  we  shall  outline  the  starting 
equations  and  geometry,  and  then  derive  the  equations  relevant  to  the  initial 
value  problem  for  infinitestinal  perturbations. 

The  cold-fluid  equations  (with  pressure  0)  describing  the  nonrelativistic 


flow  of  a  non-neutral  pure  electron  plasma  are 

dtn  -|-  V  •  ( nv )  =  0,  (la) 

(dt-v-V)v  +  £  +  vxft  —  Q,  (lb) 

Z  =  -V<£,  (lc) 

V20  =  (Id) 


Here,  n,  m,  and  v  denote  the  electron  number  density,  the  electron  mass, 
and  the  velocity,  respectively.  The  plasma  frequency  is  u =  4z e2/m.  The 
equilibrium  we  consider  corresponds  to  the  planar  configuration  of  Fig.  1  with 
crossed  electric  and  magnetic  fields.  The  cathode  is  at  y  =  0,  and  the  anode 
at  y  —  We  define  the  normalized  electric  field  £o  =  eE0/m  =  -eEoy/m 
and  the  gyrofrequency  Q.  =  —e(Bo/rnc)z. 

In  this  paper  we  consider  electrostatic  modes  where  the  magnetic  field 
remains  equal  to  the  equilibrium  or  zeroth-order  value  at  all  orders.  Hence 
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we  will  omit  the  subscript  0  on  the  ft  even  though  this  is  a  zeroth-order 
quantity.  We  shall  also  do  the  same  for  the  plasma  frequency  shortly.  We 
consider  a  two-dimensional  diode  structure,  with  translational  invariance  in 
the  z  direction  both  in  the  equilibrium  and  for  the  perturbation  quantities. 
In  addition,  we  assume  translational  invariance  in  the  x  direction  for  the 
equilibrium  quantities.  From  Poisson’s  equation,  we  have 

dy£o  =  w2  (2a) 

where  now  u;2  =  47rn0e2/m,  and  no  is  the  equilibrium  electron  density.  The 
electron  density  n0  may  be  specified  to  be  of  any  form  and,  hence,  is  arbitrary. 
For  no  or  to2  a  monotonic  decreasing  function  of  y,  the  linear  diocotron  mode 
is  stable  for  a  cold  massless  low-density  electron  plasma.16  This  criterion  has 
recently  been  generalized  to  warm  plasmas  without  assuming  the  guiding- 
center  approximation.17  The  equilibrium  electron  flow  velocity  is 

(26) 

Equations  (2a)  and  (2b)  imply  that  u0  increases  monotonically  with  y, 
and  thus  we  have  a  nonzero  velocity  shear  in  the  equilibrium. 

We  next  consider  first-order  perturbations  of  Eqs.  (1).  In  general,  all 
physical  quantities  x  may  be  written  as 


X  =  Xo  +  £Xi  +  c2X2  H - >  (3a) 

where  e  is  a  measure  of  the  (small)  deviation  from  the  equilibrium  value  xo- 
Here,  x  may  denote  the  density,  electric  field,  or  velocity.  In  first  order,  all 
quantities  are  expanded  as 
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(36) 


Xi  =  Xi  (y,t)e,kx +  c.c. 


In  this  case,  Eqs.  (1)  reduce  to 


Dhi  4-  ikn0vix  +  dy{noV\y)  =  0 


Dv\x  —  A7vly/n  —  ik<j>  i  =  0 


(4a) 

(46) 


where 


Dviy  +  Uvix  —  dy4>i  =  0 


u>; 


{dy  ~  k  )<P\  -  ni  I  -E  I  =  0 
\n0. 


(4c) 

(4rf) 


D  =  dt  +  ikvo 


(5) 


A2  =  ft2  -  u,\ 


(6) 


One  can  simplify  the  final  equations  by  using  Lagrangian  displacements. 
If  we  define 


Uli  =  -  -jffv 


(7a) 


Uiv  = 


(76) 


fii  =  -  dv(no£v) 


(7c) 
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then  (4a)  is  identically  satisfied  and  the  remaining  three  equations  become 


DD^X  —  f ID(V  —  ik4>\  =  0 

(8a) 

DD(y  +  ClD(x  -  wl(y  -  dv4>  1  =  0 

(86) 

(dy  -  k7)4> i  +  iku}j(x  +  dy{u>l(y)  =  0 

(8c) 

The  smooth  bore  boundary  conditions  are 

-w 

H 

•O 

II 

O 

II 

cT 

r* 

(9) 

which  arises  from  the  vanishing  of  the  parallel  electric  field  at  the  cathode 

and  the  anode. 

Given  the  value  of &(y,  <  =  0),  &(y,  t  =  0),  vlx(y,  t  =  0),  and  v\ v(y,  t  =  0), 

we  wish  to  construct  the  solution  of  (8)  for  all  t  >  0.  This 
problem  we  shall  now  solve. 

is  the  initial  value 

If  we  define  the  Laplace  transform  of  a  variable  as 

Ex(y,s)=/  e~gt(x(y,t)dt 

Jo 

(10) 

by  a  capital  symbol,  then  eqs  (8)  becomes 

p2Zx  -  tlpEy  -  ik$  =  Fx 

(11a) 

P2Z„  +  ttpZ,  -  U^Zy  -  3yt  =  Fy 

(116) 

(a]  -  k2)t  +  ikwl zc  +  a,(u,7rZy)  =  0 

(11c) 
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where  $  is  the  Laplace  transform  of  4> i , 

p  =  s  +  ikvo(y) 

(12) 

and  the  inhomogeneous  terms  in  (11)  are 

Fr  =  pUt  =  0)  +  x>u(t  =  0)  -  ~iv(t  =  0) 

(13a) 

F„  =  piy{t  =  0)  +  viy(t  =  0)  +  Ct^(t  =  0) 

(136) 

By  eliminating  $  from  (11),  one  can  reduce  (11)  to  the  two  equations 


E.)  -  ik  (l  -  2^)  =,  =  ^(V  ■  F )  -  i(V  x  F),  (14a) 

+  iE,  =  -j(V  •  F)  +  ^(V  x  F),  (146) 

where 


a  =  p2  +  n2 

(15) 

(v  •  F)  =  »*fx  +  avF„ 

(16a) 

(v  x  f),  =  u-fv  -  av& 

(166) 

One  can  now  solve  the  initial  problem  by  solving  (14)  for  Ex(y,$,k)  and 
E v(y,s,fc)  subject  to  the  boundary  conditions  (9).  Through  (11a),  these 
boundary  conditions  will  relate  Ex  and  Ev  at  y  =  0  and  at  y  =  £.  This 
solution  is 
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(P2^*  —  flpHy  —  /’*)|y=0,/  =  0  (17) 

Given  the  initial  data,  the  solutions  of  (14),  which  are  inhomogeneous  equa¬ 
tions,  will  be  unique  except  possibly  at  certain  discrete  values  of  s,  which 
are  the  eigenvalues.  At  these  eigenvalues  of  s,  the  function  Ex  and  E„  satisfy 
the  homogeneous  part  of  eqs  (14).  In  actuality,  the  solution  of  the  initial 
value  problem  will  give  Ex  and  Ev  a s  having  a  simple  pole  in  s  (for  all  y)  as 
s  approaches  any  one  of  these  eigenvalues14,15. 

But  the  solution  of  (14)  also  has  other  singularities  in  s.  These  are  y 
dependent  singularities  and  occur  at  the  singular  points  of  the  second-order 
differential  system  (14).  From  (14),  it  is  obvious  that  the  points  p  =  0  and 
A  =  0  are  ordinary  singular  points  for  the  differential  system  (14).  The 
position  of  these  singularities  depend  on  s  through  (12)  and  (15).  Since  we 
shall  have  need  of  them  later,  we  shall  give  the  form  of  the  solution  near 
these  singularities. 

Near  p  =  0,  and  as  p  — >  0,  except  for  an  overall  multiplicative  constant 
(in  y,  but  possibly  dependent  on  s  and  k)  the  various  Laplace  transforms 
will  approach 


-  _  *'K  , 

“r  n  v  ' 

(18a) 

iky* 

^  Sip  +"'  . 

oo 

„  (V?) ,  , 

Vx  jp-  In  p  +  •  •  • 

(18c) 

ikul 

V»-Tf  +  - 

(1M) 
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tk 

ik$  — - - -  +  •  •  • 

(18c) 

dv$  -►  -(5yu;J)lnp  +  ••• 

(18/) 

d..u? 

N  — ♦  -ikno—r-Z-  +  •  •  • 

Up 

(185) 

where  Vx(Vy,  N)  is  the  Laplace  transform  of  0u(0iv,  Hi). 
Similarly  near  A  — »  0,  we  have 


2fi2 

(19a) 

Ev-^  +  ... 

v  A 

(196) 

+  ... 

/I 

(19c) 

i/  2n3 

+- 

(1M) 

ik$  - -tf In /!  +  •■• 

(19c) 

2pf!w* 

«/„$  — ♦ - H - 

_ 

(19/) 

w2f}2 

•N  -»  -iikno-Z-r-  +  •  •  • 

(195) 
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Now,  every  line  in  eqs  (18)  and  (19)  have  also  logarithmic  terms  in  higher 
order,  but  near  the  singularity,  they  only  dominate  in  (18c)  and  (19e).  Nev¬ 
ertheless,  their  existance  means  that  in  addition  to  the  pole  singularities 
discussed  above,  the  solutions  of  (14)  must  also  have  branch  cuts  in  the  com¬ 
plex  s-plane.  These  cuts  will  be  located  at  all  values  of  s  for  which  either  p 
or  A  will  vanish  for  some  value  of  y  between  the  cathode  and  anode.  From 
(12)  and  (15),  one  can  see  that  these  branch  cuts  would  lie  along  the  imagi¬ 
nary  s-axis,  and  in  general  would  have  three  sections  as  shown  in  Fig  2.  The 
quantity  vj  is  the  value  of  vo(y)  at  the  anode.  If  kvj  >  fl,  then  the  three 
sections  would  merge  into  one  branch  cut. 

Now  that  we  have  described  the  Laplace  transform  of  the  general  solution, 
let  us  next  look  at  the  asymptotic  form  of  the  general  solution. 

3  The  Large  t  Limit  of  the  Initial  Value 
Problem 

We  shall  now  obtain  the  large  time  asymptotic  solution  of  the  initial  value 
problem.  Given  the  Laplace  Transform,  say  Ex(y,s,k),  then  the  original 
function,  (x(y,t,k)  is  constructed  from 

i :(y  *  *1  =  2 ~i  J  ioo  5’  (20) 

where  the  contour  is  to  be  taken  to  the  right  of  all  singularities  in  the  complex 
s-plane.  The  functions  are  all  analytic  for  large  s  and  also  vanish  as  |  s  |-+  oo. 
Therefore  we  may  close  the  contour  to  the  left  of  the  imaginary  axis,  and 
shrink  it  down  into  individual  contours  around  all  branch  points  and  other 
singularities.  These  final  contours  are  indicated  in  Fig  2. 

There  are  three  basic  contributions  to  the  solution  of  the  initial  value 
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problem:  i)  contours  around  the  eigenvalues,  ii)  contours  around  the  branch 
cuts,  and  iii)  inside  each  branch  cut,  for  any  fixed  value  of  y,  there  will  be  a 
regular  singular  point  of  the  differential  system  (14).  These  latter  singulari¬ 
ties  are  tabulated  in  (18)  and  (19). 

The  contribution  from  the  eigenvalues  will  be  a  (global)  eigenmode  with 
a  frequency  of  /m(sj)  and  a  growth  rate  of  .fte(sj)  where  Sj  is  the  eigenvalue. 
This  contribution  has  been  very  well  described  in  the  literature1-8.  The 
contributions  from  the  branch  cuts  (excluding  all  pole  singularities  inside 
the  cut)  will  always  decay  at  least  algebraically  in  time9-11,14'15,  usually  at 
least  as  fast  as  t~l.  Thus  it  is  only  a  transient  part  of  the  solution. 

On  the  other  hand,  the  pole  singularities  inside  the  branch  cuts  can  give 
a  larger  contribution.  For  example,  (18a)  has  a  pole  or  order  2  about  p  =  0. 
Inserting  such  a  pole  into  (20)  will  give  a  growing  contribution.  Let  a(s,y) 
be  analytic  in  s  near  p  —  0.  Then  one  may  easily  show'  that  for  large  times 

5 <»> 

Similarly 

2~^  f  Q(s,  3/)e,t(ln  p)ds  — >  0  (23) 

From  (21  )-(23),  one  can  easily  read  off  the  asymptotic  values  of  the  various 
quantities  in  (18)  and  (19).  For  (18),  near  p  =  0,  we  see  that  for  large  t 

l\x  — ►  ikCotj£e~,kvoi  H -  (24a) 
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(246) 


£iv  —*  ikCot-^-e~ikvoi  +  •  •  • 
vix,  Civ,  ik4>\ » “*  0  (24c) 

& — ao,  (^)  (Sfli) . ...  (24.) 

where  Co  is  the  overall  normalization  constant,  Co(s,k),  for  (18),  evaluated 
at  p  =  0.  This  constant  is  determined  by  the  initial  data,  and  in  general  is 
nonzero. 

As  (24)  shows,  the  solution  is  rather  simple.  For  large  t,  the  y  displace¬ 
ments  approach  a  constant  amplitude  and  are  carried  along  with  the  back¬ 
ground  flow,  v0.  (Note  that  when  the  e~lkvot  term  is  combined  with  the 
elkx  term  in  (3b),  then  £ir,£iv,  and  hy  will  be  functions  of  the  combination 
(x  —  Vo*))-  Of  course,  any  displacements  in  the  y-direction  shifts  the  particle 
into  a  different  background  flow.  Whence  there  is  a  lateral  displacement  of 
the  fluid  particle  away  from  its  original  position  which  grows  linearly  in  time. 
This  lateral  displacement  is  also  directly  proportional  to  the  background  ve¬ 
locity  difference  between  the  initial  and  final  layer.  This  gives  £*  =  f£ydvv0 
and  is  exactly  the  result  found  by  (24a).  At  the  same  time,  the  density  shifts 
according  to  the  vertical  displacement  of  the  particles.  Such  motion  is  well 
known  from  the  fluid  dynamics  of  shear  flows18.  Nothing  more  needs  to  be 
said  about  it  here  except  that  it  verifies  that  the  treatment  so  far  is  correct. 

Now  consider  the  solution  at  the  so-called  "magnetron  instability1 ,2’12’19”, 
which  is  where  A  ~  0,  or 

p  =  ±tfl  (25) 

As  before,  we  obtain 
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6*  —  =F  iC±Sle±inte-ikvot 

(26  a) 

6V  —  C±  fie±,n<  e~ikvot 

(26  b) 

vlx  —  C±A2e±‘nte-*'fcV0< 

(26  c) 

vlv  ->  ±iC±n7e±i(He-ikvot 

(26  d) 

ikj>i  — >  0 

(26e) 

dji  -  —C±Cluj7e±int e~ikvot 

(26/) 

—  -  ikC±u'lte±inte~ikvoi 
n0  p 

(26  g) 

where  C±(s,k )  are  the  constants  of  proportionality  for  (19). 

Now  the  motion  is  more  complicated.  However  the  interpretation  is  again 
simple.  Per  (26)  the  fluid  particles  now  undergo  cyclotron  motion  with  a  ra¬ 
dius  of  |C±ft|.  They  are  also  carried  along  with  the  background  flow.  The 
resulting  fluid  velocity  is  assymetrical20,21  (see  also  (7)).  The  electrostatic 
potential  and  the  parallel  electric  field  both  vanish,  but  the  vertical  electric 
field  does  not  vanish.  Instead  it  approaches  a  constant  amplitude  propor¬ 
tional  to  the  density,  and  is  exactly  equal  to  As  a  consequence  of  this 

vertical  electric  field,  there  is  a  constant  flux  of  particles  into  this  motion, 
with  the  number  of  particles  involved  in  the  motion  growing  linearly  in  time. 
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Note  that  the  electrostatic  potential  does  vanish,  but  that  the  fluid  veloc¬ 
ities  do  not.  As  such,  it  is  the  fluid  motion  which  is  most  significant  for  this 
solution.  The  electrical  potential  has  the  least  importance.  It  vanishes.  The 
vertical  electrical  field  seems  to  arise  as  a  result  of  only  the  cycling  cyclotron 
particles.  When  they  move  up  vertically,  the  disturbed  charge  density  cre¬ 
ates  a  more  intense  vertical  electric  field  at  that  point.  And  the  resulting 
imbalance  in  the  vertical  component  of  the  force  equation  due  to  the  shear 
flow  requires  this  vertical  electric  field. 

Another  way  to  look  at  this  solution  is  to  consider  the  effect  of  electrons 
undergoing  cyclotron  motion  in  a  shear  flow.  Even  if  they  would  be  initially 
in  phase,  due  to  the  shear  flow,  different  layers  will  go  out  of  phase,  creating 
large  vertical  density  gradients  which  generate  the  vertical  electric  field.  And 
this  time  oscillating  electric  field  will  then  pull  more  particles  into  a  cyclotron 
motion.  Note  that  the  flow  is  compressible  since  V  •  v  /  0. 

This  motion  seems  to  describe  the  recent  numerical  simulations  of  a  2  vane 
magnetron  done  by  MRC22.  In  Fig  3  we  show  a  series  from  these  simulations. 
What  one  can  observe  is  a  slow  growth  of  a  circular  structure  which  extends 
over  two  vanes.  This  structure  is  actually  composed  of  rotating  particles  and 
is  called  a  convective  cell23-25.  The  number  of  particles  involved  in  the  motion 
do  seem  to  grow  linear  in  time.  Unfortunately,  exact  numerics  and  statistics 
adequate  for  a  good  quantitative  comparison  are  currently  unavailable. 

However  simulation  results  on  a  10  vane  magnetron26  indicated  that  con¬ 
vective  cell  formation  was  less  visible  and  certain.  Although  one  could  ob¬ 
serve  an  occasional  sporatic  formation  of  vorticies  which  would  seem  to  slowly 
grow  for  some  indefinite  time,  nevertheless  they  would  eventually  collide  and 
dissipate  away.  In  any  case,  their  formation  was  never  observed  as  clearly 
and  as  certain  as  was  in  the  2  vane  simulations. 


This  is  not  inconsistant  with  these  results,  because  in  the  2  vane  case,  only 
the  fundamental  A:- vector  and  its  harmonics  could  exist  due  to  the  periodicity. 
Thus  one  had,  in  effect,  only  one  value  of  k  present.  But  in  the  10  vane 
simulations,  sidebands  of  k  (0.9k,  1.1k,  etc)  could  simultaniously  coexist 
with  the  fundamental.  One  could  approximate  this  situation  with  an  almost 
continous  spectrum  of  fc,  where  the  general  solution  would  be  an  integral 
over  k.  As  a  consequence  of  this,  there  would  be  an  additional  dephasing 
which  would  reduce  the  growth1'4,15  predicted  by  (24)  and  (26).  And  this  is 
consistant  with  the  simulations  of  the  10  vane  magnetron. 

4  The  Direct  Asymptotic  Solution 

If  the  solutions  (24)  and  (26)  are  true,  then  one  should  be  able  to  obtain 
them  directly  from  (8).  Indeed  one  can.  For  the  p  — *  0  solution,  simply  take 
the  ansaltz 


7  v'  7^0)(y,k) 

*  "  6  *n+2 

n=0  1 

(27a) 

«.  =  f;  ^  —  no M 

(27  b) 

n=0  1 

±v°'Xk) 

n=0  1 

(27c) 

U»V  =  e  L  fn+2 

nsO  1 

(27  i) 

Then  eq.  (8)  provides  recursion  relations  which  determines  all  the  func¬ 
tions  7n°\  VnK  vi°\  and  u in  terms  of  just  one  function.  Choosing  this 
function  to  be  7q°\  we  find 
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(28a) 


(o)  _  (0) 

70  —  2  70 

(28a) 

(286) 

4o) = 

(2  8c) 

(o)  20  (o) 

7i  =  a 

(29  a) 

(0)  iklo0’  a  1 

“  n j>6+’’ 

(296) 

(29c) 

(0)  _  -n2  8  (0) 

1  ~  AW 

p 

(29  d) 

and  etc.  Note  the  inverse  powers  of  in  (29a, b,d).  These  terms  indicate 
that  this  asymptotic  series  will  require  a  longer  time  for  the  higher  order 
terms  to  become  smaller  wherever  the  density  is  small.  And  the  relative 
sizes  of  these  terms  indicate  the  time  required  for  a  certain  accuracy  to  be 
achieved.  Comparing  (28a)  and  (24d)  reveals  that 


(o)  _  QCp  dvn0 
ik  no 


(30) 


Since  solving  (11)  can  give  one  the  value  of  Co,  we  see  that  the  asymptotic 
solution  (27)  can  also  be  related  to  the  initial  data. 
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Now  consider  the  magnetron  instability  where  A  — ►  0.  The  correct  ansaltz 
is 


2  —  c±iClt  -ikvat 

V  fn+ 2 

n=0  1 

(31a) 

«.  _  e*n.e-««<  £;  k)Mv) 

(316) 

n=0  1 

_  c±«ntc-ifcw)t  y' 

x  4n 

n=0  1 

(31c) 

0ly  -  e±‘n<e-‘W  £  fc) 

n=0 

(31rf) 

Again,  the  lowest  order  solutions  are 


/O  ^2 

(32a) 

(326) 

=  Tfc7o:t> 

(32c) 

o 

II 

* 

(33a) 

”!±)  =  -ir}d-  M 

(336) 

,(±)  =  ±0*  <±>  _  .(±) 

1  ^7o  n^r*70 

(33c) 
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and  etc.  As  before,  convergence  will  be  poorest  in  the  low  density  regions. 
Also  from  (32a)  and  (26g) 

7?’  =  |  (34) 

Note  that  in  (26),  C±  is  the  overall  constant  of  proportionality  for  the 
solution  (19).  Thus,  C±  can  be  found  from  the  initial  data,  and  by  (34),  7^ 
can  also  be  found  from  the  initial  data.  Thus  given  the  initial  data,  one  can 
in  principle  construct  the  time  asymptotic  solution  by  knowing  simply  three 
function  of  y  and  k;  namely  7o ^(y,^)  and  7o 

5  Summary 

We  have  presented  a  generic  algebraic  instability  in  the  planar  magnetron. 
It  will  always  be  present  except  for  those  rare  initial  data  for  which  the 
constants  C±  vanish.  It  is  furthermore  independent  of  the  density  profile. 
Thus  even  if  a  stable  equilibrium  density  profile27  existed,  this  algebraic 
instability  would  cause  the  density  profile  to  drift  away  from  equilibrium. 

This  instability  seems  to  have  been  important  in  the  2  vane  magnetron 
simulations.  However,  it  seemed  to  be  less  dramatic  when  sidebands  could 
be  present,  as  in  the  10  vane  simulations. 


Figure  Captions 

Fig.  1:  Geometry  and  the  shear  flow  in  the  planar  magnetron. 

Fig.  2:  The  complex  s*plane,  showing  the  branch  cuts  (  x  -  -  -  x  )  and  the 
location  of  possible  eigenvalues  (*)  for  the  Laplace  Transforms.  At 
each  value  of  y,  there  is  also  a  singularity  inside  each  branch  cut  which 
is  located  at  —ikvo{y)  below  the  top  of  the  branch  cut.  Thus  these 
singularities  move  as  a  function  of  y.  Vj  is  the  value  of  vo(y  =  £)■  The 
contour  C  is  discussed  in  the  text. 

Fig.  3:  The  formation  of  a  convective  cell.  (Courtesy  of  Mission  Research 
Corp.) 
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FIG.  1  - 


Geometry  and  the  shear  flow 
in  the  planar  magnetron 


is  is  not  inconsistant  with  these  results,  because  in  the  2  vane  case,  only 
tdamental  A;- vector  and  its  harmonics  could  exist  due  to  the  periodicity, 
me  had,  in  effect,  only  one  value  of  k  present.  But  in  the  10  vane 
tions,  sidebands  of  k  (0.9k,  1.1k,  etc)  could  simultaniously  coexist 
he  fundamental.  One  could  approximate  this  situation  with  an  almost 
ous  spectrum  of  k ,  where  the  general  solution  would  be  an  integral 
.  As  a  consequence  of  this,  there  would  be  an  additional  dephasing 
would  reduce  the  growth1'4,15  predicted  by  (24)  and  (26).  And  this  is 
tant  with  the  simulations  of  the  10  vane  magnetron. 

The  Direct  Asymptotic  Solution 

solutions  (24)  and  (26)  are  true,  then  one  should  be  able  to  obtain 
directly  from  (8).  Indeed  one  can.  For  the  p  — »  0  solution,  simply  take 
isaltz 


i  =  e-*».  f 

V  Zrf  in+2 


(11a) 


n=0  1 

-  _  V'  vn]{y^) 

v'*~e  1 ,  -+r" 


(276) 


(27c) 


_  ^f«r(y,i) 

UlV  -  C  Z.  yrT+2 


(27  d) 


Tien  eq.  (8)  provides  recursion  relations  which  determines  all  the  func- 
ln\Vn\vn\  &nd  terms  of  just  one  function.  Choosing  this 

tion  to  be  7q°\  we  find 
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FIG.  3  -  The  formation 
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In  1975,  one  of  the  present  authors1  showed  how  to  obtain  the  quantized  levels  of 
the  nonlinear  Schrodinger  equation  using  the  action-angle  variables  (canonical  coordinates) 
of  the  AKNS  scattering  data.  The  symplectic  form  used  to  effect  the  reduction  to  canonical 
coordinates  was  based  on  the  standard  Hamiltonian  structure  for  the  nonlinear  Schrddinger 
equation.  The  method  used  was  a  nonlinear  generalization  of  one  of  the  standard  methods 
for  the  second  quantization  of  the  electromagnetic  field.  As  presented  in  the  textbook  by 
Schiff2,  one  takes  the  classical  electromagnetic  field  and  decomposes  it  into  normal  modes 
(Fourier  components).  The  key  idea  in  this  approach  is  that  the  classical  electromagnetic 
Hamiltonian  will  decompose  into  a  sum  of  noninteracting  classical  Hamiltonians,  each  of 
which  has  just  two  degrees  of  freedom  and  is  easily  quantized  by  itself.  This  method  of 
quantization  bypasses  all  the  inherent  difficulties  of  fully  quantizing  the  system,  including 
the  factor-ordering  problem,  defining  the  quantum  field  operators  for  the  fundamental 
fields,  etc.^  It  is  fundamentally  based  on  the  symmetries  of  the  classical  system,  and 
reduces  the  problem  to  one  of  quantizing  noninteracting  particles4.  In  this  way,  the  original 
difficult  second  quantization  problem  is  reduced  to  a  simpler  set  of  noninteracting 
problems.  The  advantage  of  this  simpler  solution  is  tremendous  when  one  considers  the- 
information  that  one  can  glean  from  it.  First,  one  can  obtain  the  spacings  of  the  energy 
levels.  One  also  discovers  which  quantum  variables  will  commute,  and  which  modes  will 
have  a  panicle-like  behavior.  Of  course,  for  a  full  quantum  theory,  one  still  has  to  deal 
with  a  number  of  remaining  difficult  problems,  including  finding  a  consistent  factor- 
ordering  for  the  quantum  operators,  evaluating  matrix  elements,  etc.  Unfortunately,  the 
solution  to  this  larger  quantization  problem  may  well  be  multi-valued  .  However,  in  the 
meantime,  one  has  been  able  to  immediately  isolate  the  above-mentioned  important  features 
of  second  quantization,  and,  very  importantly,  those  quantities  which  would  have  the  same 
common  solution  for  every  possible  consistent  second  quantization.  Thus,  any  difficulty 
which  would  be  found  at  this  level  would  also  be  present  in  any  quantum  field  theory.  And 
a  study  by  this  method  can  provide  valuable  insight  into  the  structure  of  the  more  thorny 
parts  of  the  second-quantization  problem. 

The  symplectic  form  used  in  ref.  1  to  effect  the  reduction  to  canonical  coordinates 
was  based  on  the  first  Hamiltonian  structure  for  the  nonlinear  Schrddinger  equation.  In 
1978,  Magri5  showed  how  many  soliton  equations,  including  the  nonlinear  Schrodinger 
equation,  could  be  written  as  biHamiltonian  systems,  meaning  that  they  have  two  distinct, 
but  compatible,  Hamiltonian  structures.  Indeed,  his  fundamental  result  showed  that, 
subject  to  some  technical  hypotheses5,6  any  biHamiltonian  system  is  completely  integrable 


in  the  sense  that  it  has  infinitely  many  conservation  laws  in  involution  and  corresponding 
commuting  Hamiltonian  flows. 

From  the  viewpoint  of  quantum  mechanics,  the  existence  of  more  than  one 
Hamiltonian  structure  for  a  given  classical  mechanical  system  raises  the  possibility  of  there 
existing  more  than  one  quantized  version  of  this  system,  even  at  the  level  of  quantization 
considered  in  ref.  1.  The  resulting  ambiguity  in  the  quantization  procedure  raises  serious 
physical  doubts  as  to  the  mathematical  framework  of  quantization.  However,  the  main 
result  to  be  proven  here  is  that,  for  AKNS  soliton  equations7,  both  quantized  versions  are 
essentially  the  same.  We  demonstrate  that,  in  terms  of  the  respective  canonical  coordinates 
on  the  scattering  data,  the  two  Hamiltonians  have  identical  expressions,  and  hence  identical 
quantum  versions.  Indeed,  we  conjecture  that  this  phenomenon  is  true  in  general  - 
quantization  does  not  depend  on  the  underlying  Hamiltonian  structure.  (The  results  of 
Dodonov  et.  al.8,  in  which  an  ambiguity  in  the  quantization  procedure  for  certain  finite¬ 
dimensional  biHamiltonian  systems  is  supposedly  demonstrated,  are  erroneous,  since  they 
fail  to  incorporate  the  important  topological  properties  of  phase  space  properly  in  their 
picture.  Indeed,  their  ambiguity  is  just  a  version  of  the  ambiguity  inherent  in  the: 
quantization  of  two-dimensional  Hamiltonian  systems,  which  we  discuss  in  detail  below.) 
Moreover,  we  will  see  that  for  the  other  members  of  the  associated  hierarchy  of  soliton 
equations  the  only  difference  in  the  quantum  versions  is  in  the  choice  of  weighting  factor 
for  the  quantum  operators  corresponding  to  the  continuous  spectrum,  the  weight  being 
determined  by  the  classical  dispersion  relation,  and  the  replacement  of  the  bound  state 
Hamiltonians.  Thus,  the  effect  of  quantizing  different  members  of  the  soliton  hierarchy 
will  only  be  significant  for  the  bound  states/solitons. 

Our  presentation  relies  heavily  on  the  notation  and  results  in  earlier  papers  by  Kaup 
and  Newell1,9,1®  on  the  closure  of  the  squared  eigenfunctions  for  the  AKNS  scattering 
problem.  The  key  to  our  result  is  the  well-known  fact  that  the  recursion  operator,  which  is 
built  out  of  the  two  Hamiltonian  operators  for  the  system5’®  is  essentially  the  squared 
eigenfunction  operator.  Since  variations  in  the  potential  for  the  AKNS  mattering  problem 
are  expressed  in  terms  of  the  squared  eigenfunctions,  this  means  that  the  second  symplectic 
form  can  be  simply  written  down  in  explicit  form;  in  terms  of  the  scattering  data,  it  differs 
from  the  first  symplectic  form  only  by  a  weighting  factor  in  the  continuous  spectrum,  and  a 
change  in  the  discrete  components.  However,  the  corresponding  difference  in  weighting 
factors  for  the  two  Hamiltonians  exactly  cancels  out  the  weighting  factor  for  the  two 
symplectic  forms,  while  the  discrete  components  reduce  simply  to  the  quantization  of  a  two 


dimensional  Hamiltonian  system,  based  on  different  symplectic  structures.  Thus,  the  entire 
quantum  ambiguity  reduces  to  the  simple  matter  of  an  ambiguity  in  the  quantization  of  two 
dimensional  Hamiltonian  systems,  a  problem  that  is  easily  handled. 

Our  notation  is  as  follows.  Hamilton's  equations  are 

9tQa  =  JaPapH,  (1) 

where  Q={Qa}  are  the  dynamical  variables  (the  p’s  and  the  q's),  J  =  [J0^]  is  the 
Hamiltonian  operator,  which  determines  the  underlying  Hamiltonian  structure  of  the  phase 
space,  and  H  is  the  Hamiltonian  function  or  density.  For  instance,  for  a  harmonic 
oscillator,  one  would  take 

Q  =  (p) .  **  H=i(p2  +  q2). 

When  Q  is  a  function  of  a  continuous  variable,  the  sum  over  the  dummy  indices  in  (1)  is 
understood  to  include  the  appropriate  integration,  and  the  partial  derivative  is  understood  to 
be  a  functional  derivative  instead.  The  Poisson  brac<et  determined  by  such  a  Hamiltoniaiv 
operator  has  the  form 

{ F,  G  }  =  (90  F)  9jj  G  ,  (2) 

which  requires  the  symplectic  two-form  to  be 

ft  =  i  dQaAj^dQP.  (3) 

For  the  harmonic  oscillator,  this  reduces  to  the  familiar  canonical  form 

ft  =  dp  a  dq.  (4) 

Therefore,  the  operator  J  needs  to  be  skew  adjoint,  and  satisfy  the  additional  condition 
that  the  Poisson  bracket  (2)  satisfy  the  Jacobi  identity,  which  is  equivalent  to  the 
requirement  that  the  two  form  ft  be  closed®. 

Before  presenting  the  main  results,  we  discuss  a  simple  but  crucial  fact  that  any  two 
dimensional  Hamiltonian  system  has  a  unique  quantized  version,  even  though  it  has  many 
different  Hamiltonian  structures.  In  terms  of  the  standard  Hamiltonian  structure  prescribed 
by  the  canonical  two-form  (4),  Hamilton’s  equations  take  the  classical  form1 1 


9H 

p'  =  ~* ■ 


9H 

q‘ '  ¥ ' 


(5) 


In  IR2,  any  nonzero  two-form  X(p,q)  dp  a  dq  is  always  closed,  and  hence  determines  a 
Hamiltonian  operator 


r  = 


l  \ 

‘x 

o 


It  is  easy  to  see  that  (5)  can  be  written  in  Hamiltonian  form  using  this  second  Hamiltonian 
structure  if  and  only  if  X  is  a  function  of  the  Hamiltonian  H.  In  this  case,  the  new 
Hamiltonian  function  is 


H2(p,q)  =  d>[  H(p,q)  ]  , 

where  C>(4)  is  any  nonvanishing  scalar  function,  and 

n2  =  C>'[H(p,q))dp  Adq  (6)  *c 

is  the  second  symplectic  form.  Re-expressing  Q2  *n  canonical  form  will  lead  to  new 
canonical  variables  p,  q,  and  an  ostensibly  different  quantized  version.  However, 
provided  this  transformation  does  not  affect  the  phase  space  topology,  it  is  not  hard  to  see 
that  these  two  quantized  versions  will  end  up  being  identical,  at  least  in  the  semi-classical 
limit,  and  so  there  is  no  ambiguity  in  the  (semi-classical)  quantization  of  two-dimensional 
Hamiltonian  systems. 

We  now  turn  to  our  problem  at  hand.  For  simplicity,  we  will  consider  the  general 
nonlinear  Schrodinger  equation 


H  - -<Ux  +  2r<iJ: 

(7a) 

'h  *  r«  +  2<ll2- 

(7b) 

in  detail.  However,  our  arguments  will  work  equally  well  for  any  other  soliton  equation 
associated  with  the  AKNS  spectral  problem7;  see  the  remarks  at  the  end  of  the  paper.  For 
r  =  ±  q*,  (7)  reduces  to  the  single  equation 

ick  =  -q*x  ±  2(q*q)q.  (g) 
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which  is  the  form  of  the  nonlinear  Schrodinger  equation  in  which  all  physical  constants, 
e.g.  h,  m,  etc.,  have  been  set  equal  to  1.  According  to  Magri5  the  nonlinear  Schrodinger 
equation  can  be  written  as  a  biHamiltonian  system 

Yt  =  Jj  9Hj  =  J23H2 .  (9) 

The  first  Hamiltonian  can  be  identified  with  the  (signed)  energy 

H1=±E=f0O(qxrx  +  q2r2)dx,  (10) 

»  oo 

while  the  second  Hamiltonian  is  the  field  momentum 

H2  =  P  =  i  J°°  (rqx-qrx)dx .  (11) 


(12) 

r- 

(13) 

(In  our  notation6  we  have  omitted  the  delta  functions  used  by  some  authors.)  Moreover, 
these  Hamiltonian  structures  are  compatible,  in  the  sense  that  any  linear  combination 
clJj  +  c2J2  is  also  Hamiltonian.  Therefore,  according  to  the  theorem  of  Magri  the  operator 

R  =  J2  •  J]1  (14) 


The  two  Hamiltonian  operators  are  given  by 

J1  -  *  (  ?  V  )  • 
'  “JI'i  -iLr  A 


h  *  2  °iax  + 


-rL^  rllr 


-(Vo) 


j 


is  a  recursion  operator  for  the  general  nonlinear  Schrbdinger  equation,  leading  to  an  infinite 
hierarchy  of  mutually  commuting  biHamiltonian  flows. 

To  determine  the  two  quantized  versions  of  the  nonlinear  Schrodinger  equation,  we 
need  to  introduce  canonical  coordinates  and  momenta,  which  will  be  found  among  the 
scattering  data  for  the  associated  eigenvalue  problem.  We  begin  by  recalling  how  this  was 
done  in  ref.  1  for  the  first  symplectic  form.  The  general  nonlinear  Schrodinger  equation  can 
be  solved  using  the  AKNS  eigenvalue  problem 

Vl,x  +  Kvi  =  qv2.  v2^ - >  C v2  -  rv,.  (15) 
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for  Im  £  <  0.  This  serves  to  define  the  scattering  coefficients  a,  b,  S,  b,  which  also 
satisfy 

S(C)a(C)  +  5(0  b(0  -  1.  (16) 

-  /- 

The  ratio  p(£)  =  b(^)/a(^),  £  real,  serves  to  define  the  continuous  spectrum  of  the 
scattering  data  for  (4).  The  zeros  of  a(Q  in  the  upper  half  plane  correspond  to  the  bound 
states,  and  are  denoted  as  £j  =  +  i  T]j,  j  =  Finally  let  bj  denote  the  value  of  b 

at  £j,  and  let  Pj  denote  the  residue  of  p  at  the  pole  Similar  quantities  are  defined  for 
the  eigenvalues 

In  ref.  1  it  was  shown  how  to  express  the  first  symplectic  two-form  in  terms  of  the 
scattering  data  in  the  case  r  =  ±  q*.  Tracing  through  the  calculation  there  in  the  more 
general  case,  we  find  that 

=  i  f  {5qA8r}dx 

*  — oo 

=  1  r  {8logb(5)A6log(S($)a($))}d5  +  (17) 

n 

N 

-2  X  {  S  Cj  A8logbj  +  8^A8logbj  }, 
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where  the  last  sum  is  absent  if  r  =  +  q*,  since  there  are  no  bound  states.  When  r  -  ±  q  , 
then  &(5)  =  a(4)*,  and  b(£)  =  +  b(£)*.  In  this  case  one  can  choose  canonically  conjugate 
variables  by  letting 

Aj  =  4  Tty  pj  =  -4^,  p(fc)  =  -  ^  log  I  a(4)  I , 

represent  the  momenta  (p's),  and 

Bj  =  arg  bj,  qj  =  log  I  bj  I ,  q(S)  =  arg  b(£). 

the  conjugate  coordinates  (q's)  for  the  system.  The  first  Hamiltonian  functional  is  then 
expressed  as 

•  N 

H1=±E=  -J~  log(  I  a£)  I )  d^  -  )•  08) 

n  ^  j=l 

From  this  expression,  the  quantized  form  follows  directly  as  in  ref.  1 . 

For  the  second  symplectic  form,  we  first  recognize  that  by  (12),  (13)  and  ref.  7, 

J2  =  LaJj  =  La02,  (19) 

where  LA  is  the  recursion  operator  for  the  squared  eigenfunctions.  Recall  that  the  squared 
eigenfunctions  corresponding  to  (15)  are  the  functions 


'P(C.x) 


Stc,*)2' 

kv2<C.X)2, 


We  define  the  corresponding  quantities  'Fj  for  the  bound  states  similarly.  The  key 
result10  is  that  the  recursion  operator  LA,  given  in  (19),  has  the  squared  eigenfunctions  as 
eigenstates: 


LAxF  =  C'F,  LATj  =  Cj'i/j.  (20) 

Thus  we  can  compute  the  second  symplectic  form 

n2  =  \  <  5Va  I  A  o2  (LA)_1 1 8V  >. 


Now,  according  to  (B3)  of  ref.  10, 
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sv  =  I  f "  (  8p(4)  >P(4)  -  8p(5)  **'(4) )  <14  - 

71 

N 

-  2  i  X  <  6Pj  +  Pj  6^j  *j +  6Pj  +  Pj  5^j  Xj)- 
j=l 


Therefore,  using  (20), 

(LV  SV  =  —  J°°  {  5p(4)  (La)_1  ¥(5)  -  8p(£)  (La)-1  9(5)  }  d£  - 

7t  “°° 

N 

-2i  £(8Pj  (La)-1  ^  +  Pj  SCj  (LV  Xj  +  5pj  (La)-j  9j  +  Pj  6^  (LV  Xj)- 
j=l 


I  r « f 

71  J— I  4  +  i  e 


5  p(5ms) 

^-i£ 


\ 

/ 


-2  i 


N  ( 


j  =  l 


*}J 


V;  + 


fa* 


+  8 


(El 

5. 


A 


V 


^84jXj 

*j 
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where  we  have  moved  the  integral  over  the  continuous  spectrum  off  the  real  axis  to  avoid 
the  singularity  at  £  =  0.  Therefore  the  only  difference  between  the  computation  of  £2j  and 
the  new  symplectic  ^nn  ft2  the  weighting  factors  1/5  in  the  continuous  spectrum, 
and  l/£j  in  the  discrete  spectrum.  A  similar  calculation  as  was  used  to  produce  (17)  now 
gives 

~  "  f  °°  (  5  lQg(  5(4)  a(5) )  a  5  arg  b(5)  }  “  + 

7t  J-~  q 


+  ^  b(0)  b(0)  6  log  ||0j  a  8  log  jj|jy  (21) 

N 

~  2  X  {  5  log  5j  A  8  log  bj  +  8  log  £j  A  8  log  bj  } , 
j=i 


where  the  two  complex  integrals  have  combined  to  give  the  principal  value  in  the  leading 
term,  and  extra  discrete  term  comes  from  the  associated  residues  at  the  pole  £  =  0.  When 
r  =  ±  q*,  canonically  conjugate  variables  are  provided  by  the  momenta 
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A.  =  4  arg  L,  p.  =  -4  log  I  L  I.  p(S)  =  -  ~  log  I  a(£)  I, 

j  j  j  j  Jt  q 

and  the  conjugate  coordinates 

Bj  =  arg  bj,  qj  =  log  I  bj  I ,  $(£)  =  arg  b(4), 

provided  £  *  0.  In  addition,  the  point  4  =  0  appears  separately  as  the  extra  residue  term 

in  the  expression  for  Cl2,  so  this  particular  mode  survives  the  principal  value  cancellation 

in  a  new  discrete  form.  However,  there  is  no  simple  formula  for  the  relevant  canonical 
variables  there.  Also,  in  the  case  r  =  ±  q*.  this  term  vanishes  because  5(0)  =  a(0),  and 
so  this  extra  complication  does  not  arise.  All  the  other  modes  for  the  continuous  spectrum 
are  related  according  to  the  simple  reweighting 

P($)  «$$($).  (22) 

For  the  second  Hamiltonian  structure,  the  Hamiltonian  functional  giving  the 
nonlinear  Schrodinger  equation  is  the  momentum  (1 1).  According  to  the  calculations  in 
ref.  1,  it  can  be  expressed  in  terms  of  the  scattering  data  as 

N 

h2  =  p=  ~J"tloglaft)l  d$  -  4i  (23) 

*  j=> 

Comparing  with  (18),  we  see  that,  in  terms  of  the  respective  canonical  variables,  the 
continuous  spectrum  contribution  is  exactly  the  same  weighted  sum  of  the  continuous 
canonical  momentum  variable  associated  with  the  respective  symplectic  two  forms: 

Hl :  “  j  Z>2  P(4)  dx  versus  H2  :  -  f "  %  p(£)  dx  a  -  J°°  t?  p(£)  dx. 
it  n  n  J-co 


Therefore,  the  continuous  modes  have  identical  quantizations.  (The  singular  point  £  =  0 
plays  no  role  as  both  Hamiltonians  make  no  contribution  to  this  mode.)  As  for  the  bound 
states,  we  are  reduced  to  the  case  of  a  collection  of  integrable  two-dimensional  Hamiltonian 
systems  with  different  Hamiltonian  structures.  For  the  original  symplectic  form  Qj,  the 
Hamiltonian  system  corresponding  to  the  discrete  eigenvalue  £j  has  the  form 


„  ,  ,  1  dH, 

0°g  bj)t  =  -  j  ^ 


1  dH, 

=0 
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and  similarly  for  the  eigenvalues  (We  are  just  reproducing  the  classical  calculation  of 
the  evolution  of  the  discrete  scattering  data  for  soliton  equations.)  For  the  second 
symplectic  form  t^ie  Hamiltonian  system  corresponding  to  the  discrete  eigenvalue 
now  takes  the  form 


(log  bj),  =  - 


1  3H1 

2  a  log  Cj 


(log  Cj)t 


1  gHi 

2  a  log  bj 


=  0, 


and  similarly  for  the  eigenvalues  £j.  Thus,  these  two  dimensional  Hamiltonian  systems 
are  identical,  even  though  they  use  two  different  Hamiltonian  structures: 


-  2  8  a  8  log  bj  versus  -  2  S  log  £j  a  8  log  bj . 


However,  as  we  remarked  above,  we  take  as  fundamental  the  fact  that  a  two-dimensional 
Hamiltonian  system  has  a  unique  quantization,  even  though  it  has  many  different 
Hamiltonian  structures.  Therefore  the  bound  states  for  the  nonlinear  Schrodinger  equation 
also  have  identical  quantizations.  We  conclude  that  both  Hamiltonians  lead  to  the  same 
quantized  version  of  the  nonlinear  Schrodinger  equation. 

As  a  final  remark,  we  recall  that  the  other  soliton  equations  appearing  in  the  AKNS 
scheme  can  be  written  in  the  form 

(?) ,  -  «o-A)  (?)  x . 


where  ft(£)  determines  the  linear  dispersion  relation7.  These  can  all  be  written  in 
biHamiltonian  form  using  the  same  two  Hamiltonian  structures  as  above.  An  identical 
calculation,  which  we  omit  for  the  sake  of  brevity,  will  show  that  the  two  quantized 
versions  of  any  member  of  these  AKNS  hierarchies  will  lead  to  the  same  quantum  version. 
Moreover,  it  is  not  hard  to  see  that  the  only  difference  between  the  quantized  versions  of 
two  different  members  of  the  same  soliton  hierarchy  is  in  the  weighting  factor  ft(£)  for  the 
modes  corresponding  to  the  continuous  spectrum  (with  appropriate  discrete  contributions  at 
the  points  where  Q(^)  =  0)  and  replacement  of  the  discrete  Hamiltonians  by  fi(Cj)  and 
Q(£j)  respectively.  Thus  the  only  distinction  between  the  various  quantized  versions  of  a 
soliton  hierarchy  is  in  the  weighting  assigned  to  the  continuous  modes,  and  the  replacement 
of  the  Hamiltonian  governing  the  evolution  of  the  bound  states.  Finally,  we  note  that  the 
same  considerations  will  apply  to  other  soliton  equations,  such  as  the  Korteweg-deVries 
equation,  as  the  key  fact  that  the  recursion  operator  is  the  squared  eigenfunction  operator 
remains  valid. 
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1.  Introduction 

The  inverse  scattering  transform  (1ST)  (1]  has  proven  to  be  a  very  powerful 
method,  allowing  us  to  fully  solve  the  initial  value  problem  for  large  classes  of 
nonlinear  partial  differential  equations  and  systems  of  ordinary  differential 
equations.  However,  except  for  certain  special  cases  (2-5],  it  has  not  been 
possible  to  solve  initial-boundary  value  problems  for  these  same  integrable 
equations.  The  basic  reason  for  this  problem  in  treating  the  initial-BVP  is  easy  to ' 
understand:  the  equations  which  determine  the  time  evolution  of  the  scattering 
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data  involve  terms  which  appear  to  require  an  overdetermination  of  the  problem 
16,  7].> 

The  difficulties  which  arise  are  generic  [6,  7],  but  let  us  consider  a  specific  case, 
the  forced  Toda  lattice  (8-15]  (FTL).  By  the  FTL  we  mean  a  semiinfinite  Toda 
lattice  [16]  in  which  the  leftmost  particle  is  externally  forced  to  move  in  a 
specified  way,  Q0(t).  No  particular  problems  arise  in  setting  up  the  1ST  for  this 
problem  (12,  17,  18].  Unfortunately,  however,  the  differential  equation  which 
determines  the  time  evolution  of  the  scattering  data  for  this  1ST  [12]  involves  as  a 
potential  the  motion,  Qx(t),  of  the  second  particle  in  the  lattice,  which  is  not 
known  until  after  the  equations  of  motion  have  been  solved.  Although  this 
situation  is  unpleasant,  it  is  not  surprising.  In  fact,  if  we  solve  the  linear 
approximation  to  the  FTL  by  Fourier  transforms  we  find  that  the  same  thing 
occurs:  the  time  evolution  of  the  Fourier  coefficients  depends  not  only  on  Q0(t), 
which  is  given,  but  also  on  Qft).  In  the  linear  case  this  apparent  need  to 
overspecify  the  boundary  condition  is  easily  resolved,  and  Qft)  can  be  found  in 
terms  of  Q0{t).  Since  it  has  proven  fruitful  in  the  past  to  think  of  the  1ST  as  a 
form  of  Fourier  analysis  for  nonlinear  systems  [1],  it  is  worth  considering  whether 
this  problem  can  be  overcome  in  the  nonlinear  case  as  well. 

Previous  work  on  the  forced  Toda  lattice  has  suggested  that  a  careful  consider¬ 
ation  of  the  analytic  properties  of  the  scattering  data  might  serve  to  resolve  this 
apparent  need  to  overspecify  the  boundary  conditions  [12,  14,  15].  However, 
these  earlier  results  were  useful  only  for  very  short  times,  since  they  involved 
calculating  the  terms  in  a  Maclaurin  series  (in  time)  for  the  scattering  data.  In 
Section  3  we  will  show  how  to  easily  recover  these  Maclaurin  series  results. 

In  Sections  4  and  5  we  establish  an  integral  kernel  representation  for  the 
scattering  data  and  show  how  the  analyticity  requirement  serves  to  restrict  the 
kernels.  We  then  use  this  restriction  to  find  a  single  nonlinear,  nonlocal  partial 
differential  equation,  involving  no  unknown  potentials,  for  the  single  independent 
kernel.  In  this  formulation  the  only  freedom  is  in  the  choice  of  the  boundary 
condition,  Q0(t).  Thus,  we  have  eliminated  the  apparent  need  to  overspecify  the 
boundary  conditions  in  order  to  determine  the  time  evolution  of  the  scattering 
data.  The  appeal  of  this  approach  is  that  it  seems  to  provide  a  direct,  but 
nonlinear,  pathway  from  specification  of  the  boundary  condition  to  calculation 
of  the  scattering  data. 

Of  course,  one  of  the  major  advantages  of  the  1ST  approach  has  always  been 
the  fact  that  it  gives  the  solution  to  nonlinear  problems  through  linear  means. 
This  advantage  is  lost  in  the  approach  discussed  above.  In  Section  7  we  present 
an  approach  in  which  the  kernel  from  which  the  scattering  data  is  constructed  is 
recovered  from  the  solution  to  a  linear  integral  equation.  In  writing  down  this 
linear  integral  equation  we  are  free  to  choose  only  one  function  of  one  variable, 
i/t(r).  The  scattering  data  for  the  forced  Toda  lattice  is,  of  course,  a  function  of 
two  variables.  Thus  we  have  found  a  reduction  of  the  problem  from  two  variables 
to  one.  No  other  method  has  achieved  such  a  reduction.  The  remaining  unsolved 


'Fokas  (26]  has  recently  established  a  nonlinear  integral-differentia]  equation  (with  singular  kernel) 
tor  the  time  evolution  of  the  scattering  data  of  the  nonlinear  Schr5dinger  equation. 
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problem  is  to  understand  the  mapping  of  the  boundary  condition,  0„ ,  into  the 
function  f/,.2  Of  course,  that  may  be  the  most  difficult  part  of  the  problem.  If  so, 
we  have  at  least  delineated  what  the  remaining  problem  is,  and  have  given  the 
analytic  structure  of  the  solution. 


2.  Background 

The  forced  Toda  lattice  [8-15]  is  a  semiinfinite  Toda  lattice  [16]  with  the  leftmost 
particle  driven  by  some  external  forcing.  The  equations  of  motion  are 

Q"  -  exp(c»-i-e„)-exp(en-eB+1),  «  =  1,2,3 .  (1) 

(Here  and  in  what  follows  time  derivatives  will  be  denoted  by  primes.)  Q„(i)  is 
the  displacement  of  the  nth  particle  in  the  lattice.  0o(f)  is  a  specified  boundary 
condition,  and  the  initial  data  (0„(Q),  0'(O)}  are  given. 

The  scattering  data  for  the  forced  Toda  lattice  can  be  recovered  from  a 
function  x(*»0  which  satisfies  the  second  order  ordinary  differential  equation 
[12] 


x"  +  [exp(  00  -  0i )  ■ +  i  Qo  -  ( i  0o  +  *  )2]  X  -  0.  (2) 

A  is  given  in  terms  of  the  spectral  parameter  z  by 


Here  Q0(t)  is  the  position  of  the  leftmost  particle  in  the  lattice,  which  is  the 
externally  specified  boundary  condition.  0,(0  is  the  position  of  the  second 
particle  in  the  lattice,  which  is  of  course  unknown  until  after  the  problem  is 
solved.  Two  different  approaches  have  been  used  to  try  and  resolve  this  diffi¬ 
culty.  Numerical  studies  [13]  have  been  used  to  make  physically  reasonable 
approximations  to  the  function  0,(f).  This  approximation  is  then  used  in 
Equation  (2)  to  find  x(*.  O-  The  other  approach,  which  we  will  use  here,  exploits 
the  analytic  properties  of  the  inverse  scattering  transform  for  the  forced  Toda 
lattice  [12,  14, 19,  20], 

Analysis  of  the  forced  Toda  lattice  1ST  shows  that  the  function  /?(;,/) 
defined  by 


*(*>')  -  7X(2,')exp[- i(r- j)r]  (4) 

must  be  an  analytic  function  of  z  inside  the  unit  circle,  |r|  -1  [12, 15].  In  general. 


1The  inverse  of  this  mapping  is  trivial.  Given  H,.  we  can  determine  the  boundary  condition  and  the 
entire  solution  by  solving  our  linear  integral  equation. 
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the  solutions  of  Equation  (2)  will  not  have  this  property.  Hansen  and  Kaup  (14] 
were  able  to  show  that  if  the  solutions  x(*.  0  Equation  (2)  are  required  to 
have  this  analytic  property,  then  the  Maclaurin  series  (in  t)  of  the  unknown 
potential  Q^t),  and  hence  the  Maclaurin  series  of  the  scattering  data  x(*.  0.  are 
determined.  This  result  is  clearly  of  only  limited  practical  use,  since  it  only  allows 
us  to  approximate  the  scattering  data  for  very  short  times.  In  what  follows  we 


will  show  how  to  make  progress  toward  a  more  useful  result. 

Let  us  rescale  the  quantities  in  Equation  (2)  in  the  following  way: 

i  =  ~2  iz,  (5a) 

S(S,t)  -JX(M).  (5b) 

q{t)  =  -26(0.  (5°) 

r(t)  -  exp[eo(0  -  2,(01 -1  +  :2o(0  -  i[26(0]2-  (5d) 

Equation  (2)  then  becomes 

S"+  [£J  +  /*?(r)  +  r(r)]S  «  0,  (6) 

where  E  and  k  are  given  by 


Finally,  we  note  that  the  function  R(z,t)  appearing  in  Equation  (4)  becomes 

*(*,0  -  S(S,t)e~'E'.  (8) 

Therefore,  S(H,t)e~'E'  must  be  an  analytic  function  of  f  for  |f |  <  2.  The 
potential,  q(t)  appearing  in  Equation  (6)  is  given  in  terms  of  the  known  external 
forcing  2o(0  by  Equation  (5c).  r(t)  is  unknown ,  since  it  involves  Qt(t),  as  shown 
in  Equation  (5d). 

Equation  (6)  has  been  extensively  studied  for  re  (-00,00)  and  with  both 
potentials  q(t)  and  r(t)  given  (21-24J.  Here  we  will  study  Equation  (6)  for 
t  e  [0, 00)  with  initial  data 


S(£,0)  -1, 

S'(f,0)  =  iE  +  MO). 


(9) 


The  initial  data  (9)  correspond  to  a  lattice  which  is  initially  static,  except  for  the 
zeroth  particle,  which  is  impulsed  with  velocity  -  ^(0). 
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3.  Madaurin  series  for  the  scattering  data 

Hansen  and  Kaup  [14]  have  shown  that  the  requirement  that  the  scattering  data 
have  the  correct  analytic  properties  is  sufficient  to  determine  a  Maclaurin  series 
expansion  (in  time)  for  the  scattering  data.  Let  us  see  how  this  result  follows  very 
simply  from  the  initial  value  problem  (6),  (9)  and  the  analyticity  requirement. 
Define 


/•(£. 0  “  S{S,t)e-iE-Q{'\  (10) 

where 

Q(t)  -  jjf?(*)<&- 

F  is  an  analytic  function  of  {  for  ]£|  <  2,  so  we  can  write 

ftt.o-  E  *»(')( $)"•  in  <2.  (12) 

n-0  '  ' 

Using  this  expansion  in  the  initial  value  problem  (6),  (9)  we  find 
Fo'(t)  -  0, 

Fx’(t)  ~  R'0)F0(t)  +  q(t)F0'(t)  +  F0"(t), 

f;  =  R’F^X  +  qF;.x  +  f;lx  +  f;.2  +  qFn. „  nz  2, 

with  initial  data 

^n(O)  - 

f;(o)  ■=  o, 

where  /?(/)  is  given  by 

R(t)  -  r(r(J)  +  i‘?'(J)'t'iI?(0]2)^-  (15) 

•'o 


(14a) 

(14b) 


(13a) 

(13b) 

(13c) 


The  differential  equations  (13)  together  with  the  initial  data  (14a)  alone  would 
determine  F„(t)  for  n  -  0,1,2, 3,...  if  *(/)  were  known.  Since  R  contains  r,  it  is 
in  fact  unknown,  and  the  second  initial  datum  (14b)  serves  to  determine  the 
Maclaurin  series  for  /?(/). 

Solving  (13a)  subject  to  the  initial  data  (14a),  wc  find 


*o(0  “I- 


(16) 
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The  initial  data  (14b)  are  then  automatically  satisfied.  Using  (16)  in  (13b), 
together  with  the  initial  data  (14a),  we  find 

Fi(t)  =  *(')•  (H) 

Requiring  that  (14b)  also  hold  for  n  =1  then  gives 

R\ 0)  =  0.  (18) 

Similarly,  solving  (13c)  for  n  =  2,  subject  to  initial  data  (14a),  gives 

Flit)  =  f'[R'(s)R(s)  +  q{s)R’{s)+R"{s)  +  q(s)}ds  (19) 

Jo 

Requiring  (14b)  to  hold  for  n  =  2  then  gives 

R"(0) +  q( 0)  =  0.  (20) 

Continuing  this  process  will  generate  all  the  coefficients  R,n>( 0)  in  the  Maclaurin 
series  for  the  unknown  potential  R(t),  and  thus  generate  the  Maclaurin  series  for 
the  scattering  data.  To  proceed  beyond  the  Maclaurin  series  approach  we  must 
exhibit  in  a  more  explicit  way  the  dependence  of  the  scattering  data  on  f. 


4.  Integral  kernel  representation  for  the  scattering  data 

S(f,  r)  satisfying  Equation  (6)  with  initial  data  (9)  can  be  expressed  in  terms  of  a 
pair  of  transformation  kernels  N(t,s)  and  L(t,s)  as  [15] 


=  e0(,)|e,£'  +  J'  |Af(/,s)+ yL(r,s)Je 

(21) 

The  kernels  N  and  L  satisfy  a  system  of  linear  partial  differential  equations 

+  9(0(3,  +  d,)+R'{t)jL(t,s)  +  q(t)N(t'S)  =  0, 

(22a) 

[a? -a; 

+  q(t)(d,~  d,)+  R'(t)]N(t,s)  +  q(t)L(l,s)  =  0 

(22b) 

with  boundary  conditions 

*(m)  -  o, 

(23a) 

L(t,t)  -  0, 

(23b) 

1 

1) 

1 

(23c) 

(23d) 
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If  both  q(t)  and  r(t)  were  given,  the  partial  differential  equations  (22)  with 
boundary  conditions  (23)  would  uniquely  determine  N  and  L.  In  our  case  q(t)  is 
given  but  r(r),  which  appears  in  both  the  partial  differential  equations  and  the 
boundary  conditions,  is  unknown.  We  will  show  below  that  by  requiring  S(f,  t ) 
to  have  the  correct  analytic  properties  we  will  obtain  a  relationship  between  the 
kernels  N  and  L  which  serves  to  eliminate  the  unknown  quantities. 


5.  Analyticity  requirement 

As  discussed  above,  S(£,  t)e~‘El  must  be  an  analytic  function  of  £  inside  the 
circle  |?|  =  2.  In  general,  the  solutions  S(f,  t)  constructed  from  Equations  (22) 
and  (23)  via  equation  (21)  will  not  have  this  property.  Therefore,  we  require 

<6  rS(Lt)e~i£,d!  -  0,  n  «  0,1,2 .  (24) 

Using  the  representation  (21)  of  S(£,  t)  in  Equation  (24)  and  (25] 

"  4Wy"+l(/?)’  (25) 

where  J„  is  a  Bessel  function  of  the  first  kind,  we  find 

('  [N(t,s)Jn(t  +  s)  + L(t,s)Jn  +  y(t  +  s)]  ds  «=  0,  n  =  0,1,2 .  (26) 

J  -  t 

Thus  N  and  L  are  not  independent.  Using  [25] 

^  foJAu~x)Ji(x)^-,  (27) 

we  find  that  if  S(f, /)  is  to  have  the  proper  analytic  structure  in  f,  N  must  be 
given  in  terms  of  L  by 


-  -  J'L{t,u)^~-~du. 


(28) 


Equation  (21)  then  gives  the  scattering  data  in  terms  of  L  alone: 

Stt.O  -  eO<'>{e,£'  +  /' JK(r  +  r,E)+f  ]/.(/,*)*-**},  (29) 
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where 


c(2,£)  =  -  j'J-^-e‘Exdx. 

Jq  a 


(30) 


6.  A  partial  differential  equation  for  the  scattering  data 
involving  no  unknown  potentials 

We  can  use  the  relationship  (28)  between  N  and  L  in  the  partial  differential 
equations  (22)  and  the  boundary  conditions  (23)  to  find  a  p.d.e.  and  boundary 
conditions  for  L  alone  which  involve  only  the  known  potential  q(t). 

First,  notice  that  the  boundary  condition  (23a)  for  N  is  automatically  satisfied 
by  (28).  Using  (28)  in  the  boundary  condition  (23c),  we  find 

R(l)  -  2  f  (31) 

Using  (28)  and  (31)  in  the  p.d.e.  (22a)  we  find  an  equation  for  L  alone: 


+  )(d,+  d,)+  RV)]L(t,s)  -  q(t)fduL(t,u)^~r~-  =  0 


(32) 


with  the  boundary  conditions 


L(t,t)  =  0, 

L(t,-t)  =  -  |  +  \e-wn 


(33) 


Equation  (22b)  is  then  satisfied  identically,  due  to  properties  of  the  Bessel 
functions. 

Neither  the  p.d.e.  (32)  nor  the  boundary  conditions  (33)  involve  any  unknown 
potentials,  since  R(t)  is  given  in  terms  of  L  by  Equation  (31).  Thus,  given  a 
function  q(t)  (which  represents  the  external  forcing  of  the  Toda  lattice),  L(t,s ) 
and  hence  the  scattering  data  S(f,  t)  should  be  determined.  The  solutions  to  (32) 
would  thus  provide  a  mapping  from  boundary  condition  to  scattering  data. 
Unfortunately,  the  p.d.e.  (32)  is  nonlinear  and  nonlocal,  it  is  far  from  clear  how 
to  solve  it  even  for  a  simple  forcing,  and  the  existence  and  uniqueness  of 
solutions  remains  an  open  question.  To  try  and  avoid  these  problems  we  will  next 
consider  the  linear  integral  (Ge!’fand*Levitan-Marchenko)  equations  satisfied  by 
N  and  L. 
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7.  Linear  integral  equations  for  the  kernels  N(t,s)  and  L(t,s) 
S(£,  t)  satisfies  the  linear  dispersion  relation  [15] 


S(S,t)e-iEl  -  e~Q{,) 


i  r  p(r)s(r.«)«-|C' ... 


-M 


where 


Sa.t)  *  s(- y,r).  (35) 

The  contour  A  passes  above  all  the  poles  of  the  scattering  data,  p(f ).  in  the 
complex  £  plane.  (For  details,  see  reference  [15].)  Using  the  representation  (21)  of 
S(£,  f)  in  the  linear  dispersion  relation  (34),  we  find  the  following  linear  integral 
(GLM)  equations  for  N  and  L: 

f~‘[N(l,s)H2{-s-2)  +  L(t,s)H3{-s-2)]ds  ~  0, 

J  -  t 


L(t,2)~  H^t-z)-  (  '[N{t,s)Hl(-s-z)  +  L(t,s)H2(-s- *)]  ds  =  0, 

J  -  t 


where 


».<*>  ■  cr/if 


The  {#„(*)}  are  related  by 

(38) 

so  two  of  them  may  be  chosen  independently,  corresponding  to  the  fact  that  in 
the  original  ordinary  differential  equation,  (6),  there  are  two  independent  poten¬ 
tials,  q(r)  and  r(t).  In  what  follows  we  will  show  that  by  requiring  S(f,  t)  to 
have  the  correct  analytic  properties  we  reduce  to  one  the  number  of  functions 
Hn(z)  which  can  be  independently  specified. 

Using  the  relationship  (28)  between  N  and  L  in  the  GLM  equations  (36),  we 
find  that  the  {  Hn }  must  be  related  by 


H2{z)  -  ~HM, 
H3{z)  -  ~H2(z), 


(39b) 
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where 


G(z)  B  /  G(z-x)J,(jc)-j.  (40) 

Jq  X 

Thus,  once  Hx  is  given  the  other  H„  are  all  determined.  The  relationship 

•=  2H{(z)  (41) 

required  by  Equation  (38)  then  holds  identically,  due  to  properties  of  the  Bessel 
functions.  Using  Equations  (28)  and  (39a)  in  the  GLM  equation  (36b),  we  find 
(after  a  change  in  the  order  of  integration)  a  single  linear  integral  equation  for 
L(l,z): 


L(t,z)  -  //,(/  -  z)  +  fr(/-:,j  +  z)L(t,s)  ds  -  0,  (42) 

t 

where  the  kernel  T  is  given  in  terms  of  Hx  by 

r(a,p)  m  dx.  (43) 

Thus,  in  Equation  (42)  we  are  free  to  specify  only  one  function,  //,(z). 

This  procedure  for  finding  the  time  evolution  of  the  scattering  data  S(f,  f)  of 
the  forced  Toda  lattice  can  thus  be  summarized  as  follows: 

1.  Choose  a  single  function  //,(z). 

2.  Calculate  F(a,/3)  from  Hx(z)  using  Equation  (43). 

3.  Solve  the  linear  integral  equation  (42)  for  L(t,s). 

4.  Calculate  which  forcing  0(f)  led  to  this  solution  by  using  the  boundary 
condition  (23d)  on  L(t,s). 

5.  Use  Equation  (29)  to  construct  the  scattering  data  S(f, /)  from  L(t,s) 
and  0(f). 

Unfortunately,  we  do  not  know  until  after  the  fact  which  forcing  led  to  our 
solution,  since  0(f)  [and  hence  ^(f)]  is  not  known  until  after  the  linear  integral 
equation  (42)  for  L(r,  z)  is  solved.  However,  as  shown  in  an  appendix,  if  we  are 
given  a  forcing  0(f),  we  can  calculate  L  iteratively  from  the  linear  integral 
equation  (42)  and  the  boundary  condition  (23d). 


8.  Summary  and  conclusions 

We  have  outlined  three  approaches  to  the  problem  of  finding  the  time  evolution 
of  the  scattering  data  for  the  forced  Toda  lattice.  All  three  approaches  exploit  the 
analytic  properties  of  the  Toda  lattice  inverse  scattering  transform. 
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The  first  approach,  using  Maciaurin  series,  has  the  obvious  advantage  of 
simplicity,  since  it  deals  directly  with  the  ordinary  differential  equation  satisfied 
by  the  scattering  data.  However,  the  results  obtained  by  this  method  are  of 
practical  use  only  for  very  short  times. 

The  second  approach,  as  outlined  in  Section  6,  gives  us  a  single  partial 
differential  equation  to  solve,  with  the  external  forcing  displayed  explicitly.  This 
indicates  that  specification  of  the  boundary  condition  is  sufficient  to  determine 
the  scattering  data  for  the  forced  Toda  lattice.  However,  this  p.d.e.  is  so 
unpleasant  mathematically  tha;  it  is  hard  to  imagine  that  it  gives  any  practical 
advantage  over  using  direct  numerical  integration  of  the  Toda  lattice  equations  of 
motion. 

The  third  approach,  outlined  in  Section  7,  involves  the  solution  of  a  linear 
integral  equation  and  so  is  probably  the  method  which  would  be  most  useful  in 
practice.  However,  it  has  the  disadvantage  that  it  appears  to  be  impossible  to 
determine  which  external  forcing  is  being  used  until  after  the  linear  integral 
equation  has  been  solved.  This  problem  is  avoided  in  the  approach  mentioned 
above,  which  deals  directly  with  the  partial  differential  equations  (22).  However, 
the  approach  outlined  in  Section  7  has  the  advantage  of  allowing  the  scattering 
data  to  be  constructed  through  the  solution  of  linear  integral  equations,  while  the 
approach  using  the  partial  differential  equations  is  nonlinear.  In  our  view,  the 
outstanding  problem  is  to  determine  the  mapping  from  the  boundary  condition 
Q0  to  kernel  which  appears  in  the  linear  integral  equation  (42). 


Appendix.  Iterative  solution  of  the  linear  integral  equation  for  the 

scattering  data 

The  linear  integral  equation  (42)  can  be  solved  by  iteration  in  powers  of  the 
external  forcing  Q(t).  In  this  appendix  we  present  the  results  of  such  an 
expansion  through  second  order.  We  expand  L  and  its  boundary  conditions  as 


L(/,s) 


E  Lln)(t,s), 


n  —  0 


-  0, 


1°’ 

1  (-2Q(t)]" 

1  2  n\ 


n  =  0, 

ft  *1,2,3 . 


(Al) 


We  also  expand  Hl  (and  thus  also  T): 


e  /fn*). 


B-0 


(A2) 
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In  zeroth  order  we  find 


=  0, 

L(0)(t,z)  *  0. 

First  order  gives 


(A3) 


(A4) 


Finally,  second  order  gives 


Hf 2)(r)  -  4 </> Sq/2<1w "2^T~2^ ~ Q ( w ) ■ • 


L<J>(/,r)  -  //<»(  t  -  i )  -  4 /'du /(<  ~  />/2dw g ( « )  2^^-  <?<»>• 


(A5) 


'o  •'o 


The  iteration  can  clearly  be  continued  to  any  order  in  the  external  forcing  £?(/). 
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The  Third-Order  Singular-Perturbation  Expansion 
of  the  Planar  Cold-Fluid  Magnetron  Equations 

By  D.  J.  Kaup  and  Gary  E.  Thomas 


The  singular-perturbation  expansion  of  the  plasma  cold-fluid  equations  for 
crossed  fields  in  a  planar  geometry  is  presented.  The  general  expansion  is  carried 
out  to  third  order.  Various  instabilities  that  occur  in  the  first,  second,  and  third 
orders  are  discussed. 


I.  Introduction 

Magnetrons  and  other  crossed-field  devices  have  been  in  existence  since  the 
1940s  [1].  Over  the  decades  a  broad  technology  base  for  building  and  designing 
such  devices  {2J  has  been  developed.  However,  most  of  this  technology  is  a  result 
of  employing  empirical  methods.  The  major  reason  for  this  is  that  the  strongly 
nonlinear  nature  of  the  basic  crossed-field  interaction  mechanism  has  made  it 
very  difficult  to  obtain  a  good  understanding  of  the  physical  processes  which 
occur  in  these  devices.  Without  this  understanding,  accurately  predicting  the 
operation  of  crossed-field  devices  has  not  been  possible. 

Motivated  by  the  need  to  predict  and  improve  the  basic  performance  of  these 
devices  for  radar  and  other  applications,  a  novel  approach  to  developing  a 
strongly  nonlinear  theory  was  undertaken.  This  approach  was  to  apply  soliton 
techniques  to  modeling  the  nonlinear  interaction  mechanism  in  crossed-field 
devices  (3-5].  Given  that  it  would  be  extremely  difficult  to  solve  the  fully 
nonlinear  problem,  a  model  problem  was  chosen  which  could  contain  enough 
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truth  to  allow  some  predictions  to  be  made  that  could  be  checked  against 
available  experimental  data.  The  early  work  [3-5]  showed  close  agreement 
between  predictions  of  the  strongly  nonlinear  theory  and  available  data,  implying 
that  the  basic  approach  was  indeed  valid.  As  a  consequence,  extensive  effort  has 
been  undertaken  to  place  this  novel  approach  on  a  rigorous  footing  by  working 
to  bridge  the  gap  between  linear  theory  and  this  new,  strongly  nonlinear  theory 
through  the  use  of  standard  multiscale  perturbation  methods,  in  addition,  it  was 
hoped  that  assumptions  used  in  the  original  work  could  either  be  further  justified 
or  removed  entirely. 

In  this  paper  the  geometry  used  for  the  multiscale  perturbation  expansion  is 
the  planar  magnetron  shown  in  Figure  1.  This  device  has  a  cathode  at  y  *  0  and 
a  parallel  anode  aty*/,  with  a  large  positive  voltage  difference  V  applied  across 
the  anode-cathode  interelectrode  region.  In  addition,  one  has  a  perpendicular 
magnetic  field  which  is  sufficiently  large  to  return  any  emitted  electrons  to  the 
cathode,  thereby  producing  “magnetic  insulation.”  We  shall  ignore  any  z-depen- 
dence,  and  assume  that  the  device  is  infinite  in  extent  in  the  z-direction  and  that 
the  solution  is  invariant  under  any  z-translations. 

In  actual  devices,  there  is  also  either  a  “slow-wave  structure”  or  a  set  of  vanes 
on  the  anode.  The  purpose  of  these  structures  is  simply  to  select  and/or  limit  the 
possible  wavenumbers  which  will  propagate.  For  simplicity,  we  shall  assume  the 
anode  to  be  smooth,  realizing  that  any  solution  we  find  could  be  either  limited  or 
enhanced  depending  on  these  structures. 


A 

z 


Figure  1.  Geometry  and  the  shear  (low  in  the  planar  magnetron. 
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At  startup,  one  first  turns  on  the  magnetic  field,  then  the  voltage.  What  is 
thought  to  occur  is  a  rapid  buildup  of  an  electron  sheath  around  the  cathode.  The 
final  sheath  is  considered  to  be  a  laminar  shear  flow  of  electrons  caused  by  the 
ExB  drift.  One  expects  the  number  density  of  the  sheath  to  increase  until  the 
electric  field  vanishes  at  the  cathode  due  to  the  space  charge  above  the  cathode. 
When  this  occurs,  then  no  more  current  can  be  drawn  from  the  cathode,  and 
presumably  an  equilibrium  is  established.  However,  at  first  there  is  nothing  which 
specifies  or  determines  the  density  profile  of  this  sheath.  It  is  usually  assumed 
that  the  density  profile  is  box-shaped,  as  shown  in  Figure  2(a),  and  that  the 
density  is  a  maximum  (u*  **  iij).  However,  there  is  nothing  in  the  zeroth-order  or 
the  first-order  cold-fluid  equations  to  prevent  other  density  profiles  from  occur¬ 
ring,  except  that  density  profiles  with  positive  density  gradients  are  known  to  be 
strongly  unstable  [6-8]  and  tend  to  relax  to  profiles  without  positive  density 
gradients.  Thus  the  combination  box-and-ramp  profile  shown  in  Figure  2(b)  is 
also  a  valid  profile  to  consider. 

After  the  sheath  is  formed,  what  happens  next  is  not  entirely  clear.  The 
classical  explanation  is  that  a  high-wavenumber,  linear  instability  exists  [9].  Due 
to  thermal  noise,  this  instability  will  be  excited,  will  grow,  and  eventually  will 
saturate.  While  there  seems  to  be  nothing  wrong  with  this  explanation,  neverthe¬ 
less  nobody  has  yet  built  a  device  based  on  it.  Neither  does  this  explanation 
explain  the  narrow  operating  range  of  such  devices.  This  indicates  that  either  the 
theory  is  wrong,  or  essential  physics  has  been  omitted. 

In  this  paper,  we  describe  the  singular  perturbation  expansion  of  the  cold-fluid 
equations  for  the  planar  magnetron.  Each  order  in  the  expansion  (out  to  third 
order),  is  reviewed:  we  describe  the  main  features  of  that  order,  how  it  could 
interrelate  with  the  other  orders,  and  what  its  consequences  are,  and  we  list  any 
important  open  questions.  The  purpose  of  this  is  simply  to  detail  what  is  certain, 
to  clarify  what  is  unknown,  and  to  attempt  to  specify  what  must  be  determined. 


(a)  (b) 


Figure  2.  Two  standard  density  profiles:  (a)  Brillouin  flow  (square  profile)  and  (b)  square-plus-ramp 
profile. 


60 


D.  J.  Kaup  and  Gary  E.  Thomas 


We  start  with  the  cold-fluid,  nonneutral,  single-species,  plasma  equations 
coupled  to  Poisson’s  equation: 


d,n  +  V  •  (nv)  =  0,  (la) 

d,v  +  (v-v)v-  V<f>  +  v  x  G  =  0,  (lb) 

V24>  -  n  =  0,  (lc) 


where  <f>  is  (e/m)  times  the  electrostatic  potential,  n  is  the  plasma  frequency 
squared,  v  is  the  velocity,  and  Q  =  —  Sk,  where  A  is  the  electron  cyclotron 
frequency.  In  lowest  order,  the  system  shall  be  independent  of  x  and  t,  with  only 
a  ^-dependence.  On  top  of  this,  we  shall  impose  a  small  rf  plane  wave,  which  will 
drive  the  entire  system,  and  we  shall  follow  the  consequences  of  this  perturbation 
out  to  third  order.  In  addition  to  the  plane  wave,  we  shall  also  assume  a 
long-wavelength  modulation  of  the  plane  wave,  parallel  to  the  anode.  Thus  we 
expand  all  our  physical  quantities  (n,<t>,v)  as 


n  =  n0  +  ent  4-  t2n7  +  e3w3  +  •  •  ■  ,  (2) 

where  e  indicates  the  amplitude  of  the  perturbing  rf  wave.  We  expand  the 
high-order  terms  as 


rtl  -  [V(**-"'>+c.c.]  +  n[0),  (3a) 

n2  ~  +c.c.]  +  [ /i 2e'(A 41  “ ’  +  c.c.]  +  (3b) 

where  k  is  the  parallel  wavenumber  and  A  is  the  frequency.  All  the  above 
coefficients  [n0,  n2,  n(20),  etc.)  are  considered  to  be  functions  of  the 

following  slow  variables: 


X  =  E*. 

(4a) 

T  =  et. 

(4b) 

t  =  e2r. 

(4c) 

t'  “  e3/. 

(4d) 

Needless  to  say,  unless  otherwise  stated,  all  quantities  are  also  functions  of  y, 
which  is  assumed  to  be  of  the  same  scale  length  as  x. 
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II.  Zeroth  order 

This  order  specifies  the  background  on  which  the  rf  wave  will  propagate.  The 
resulting  equations  are 
f  , 

v0  -  i)0i,  v0  -  dy<p0/n,  (5a) 

fyo  =  "o>  (5b) 

where  n0(y)  is  arbitrary.  Imposing  the  space-charge-limited  current  condition 
gives 


W*-0)  ■=  °*  <5c) 

which,  with  4>o(y  =  0)“0  and  the  specification  of  n0(y),  uniquely  determines 
the  solution.  Of  course,  must  be  tbe  aPPbed  anode-cathode  voltage, 

which  places  one  restriction  on  the  choice  of  n0(y).  In  general,  n0  may  also  be 
dependent  on  the  slow  variables  x.  T ,  t,  and  r'.  But  these  dependencies  will  not 
enter  until  the  higher  orders. 


III.  First-order  DC 

This  order  describes  how  the  background  responds  to  the  long  spatial  scale.  Since 
n0(y)  is  arbitrary  except  for  the  voltage  condition,  we  may  scale  n J0)  to  zero  by 
redefining  n0  +  £/i{0)  to  be  the  new  n0,  without  any  loss  of  generality.  Then  we 
have 


where  the  sources  are 


-  sg>. 

(6a) 

A2 

_  _..«>)  =  C(0) 

(6b) 

=  S[°y\ 

(6c) 

d^\0)  «=  sil\ 

(6d) 

=  —  dTnQ  —  dx(n0v0). 

(7a) 

S{?  =  -  dTv o  -  v08xv0, 

(7b) 

s{?  -  o, 

(7c) 

si?  -  o, 

(7d) 

oz D.  J.  Kaup  and  Gary  E.  Thomas 
and  wc  have  introduced 


A2  =  fi2  -  n0,  (8) 

which  is  assumed  to  be  positive  definite. 

From  Poisson’s  equation  (6d)  and  (7d),  if  the  voltage  difference  between  the 
anode  and  cathode  is  to  remain  fixed,  then  the  only  solution  is 

=  0,  (9a) 


which  with  (6c)  and  (7c)  gives 


Now  (6b)  gives  the  solution  of 


=  0. 


=  ^(dTvo+  VodxVo), 


(9b) 


(9c) 


while  (6a)  and  (7a)  may  be  integrated  and  combined  with  (9c),  which  results  in 

todTdy<l>o  +  (dy4>o)3ydx<t>o  ~  ”o^o  “  A2aC,  =  0,  (10) 

where  Ci(x,  T,  t)  is  a  constant  of  integration.  This  constant  can  be  evaluated  for 
the  fixed  voltage  condition  by  simply  integrating  (10)  from  the  anode  to  the 
cathode.  Since 


3r<f>o(>’  =  0)  ~  0  -  dT4>0{y  =  l),  (11) 


one  finds 


where 


QQ  = 


(12a) 


(12b) 


is  proportional  to  the  electrostatic  energy  stored  between  the  cathode  and  anode. 
An  interpretation  of  Cj  follows  upon  differentiating  (10).  This  gives 


dTn0  + 


a 


dxnc 


0. 


(13; 
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Obviously.  (13)  shows  that  lines  of  constant  /i„  are  conveclcd  by  the  flow 


u 


d^o  if 

Q  +J(Cl 


dX*o\ 

Q  J 


(14) 


Thus  Cj  is  a  vertical  velocity,  in  addition  to  the  ExB  drift  velocity,  which 
convects  n0,  and  arises  from  the  fixed  anode-cathode  voltage.  In  contrast  to  this, 
if  we  had  used  the  fixed  charge  conditions  on  the  anode  and  cathode,  then  we 
would  have  required  dy4>0(y  **  /)  to  be  fixed.  Assuming  n0(y  -  I)  =  0,  Equation 
(10)  would  then  give  Cj  ■=  0,  upon  evaluating  it  at  y  *=  /.  In  the  fixed-voltage  case, 
the  direction  in  which  lines  of  constant  n0  are  convected  could  be  either  up  or 
down,  depending  on  the  values  of  C,  and  dx<j>0.  In  general,  if  dxU  >  0,  the 
convection  tends  to  be  more  upward,  while  if  nearby  the  value  of  dxil  is  negative, 
then  in  that  region  the  convection  tends  to  be  more  downward.  As  a  consequence 
of  this,  we  would  expect  a  vortex  cell  to  form  wherever  this  counterflowing 
convection  occurred. 

Although  the  evolution  equation  (10)  is  fully  nonlinear,  it  still  could  be 
analyzed  in  the  weakly  nonlinear  limit  by  additional  singular-perturbation  expan¬ 
sion,  which  would  result  in  the  KdV  equation  (10-12).  A  report  on  this  is 
currently  being  prepared.  An  important  consequence  of  this  is  that  once  d  <j>0 
becomes  nonzero  and  it  is  positive  somewhere  and  negative  somewhere  else,  then 
eventually  at  least  one  KdV  soliton  will  grow  and  form. 

At  startup,  one  would  like  to  consider  the  background  to  be  smooth  and 
stationary  so  that  the  above  equation  could  be  ignored,  at  least  until  some  linear 
or  nonlinear  instabilities  have  had  time  to  grow  up  out  of  any  noise.  However, 
once  any  long-scale  variations  have  been  introduced  into  the  background,  then 
this  order  can  no  longer  be  ignored,  and  it  will  require  the  fully  nonlinear 
evolution  equation  (10)  to  describe  the  evolution  of  these  long-scale  background 
variations. 


IV.  First-order  fundamental 

This  section  covers  the  same  material  as  in  the  short-wavelength  linear  stability 
theory.  The  classic  work  of  Buneman,  Levy,  and  Lindson  [9]  is  the  usual 
reference  on  this  order.  However,  they  only  considered  the  box  profile  in  Figure 
2(a),  as  have  many  others  since  (13-15).  More  recently,  Davidson  and  others 
(6-8)  have  considered  various  other  profile  shapes  in  studies  of  the  resonant 
diocotron  instability.  However  he  has  mostly  emphasized  profiles  with  positive 
density  gradients,  since  these  demonstrate  the  strongest  resonant  diocotron 
instability. 

As  far  as  is  known  to  us,  no  linear  stability  analysis  has  ever  been  done  on  the 
ramped  profile  shown  in  Figure  2.  To  understand  why  this  latter  profile  would  be 
of  importance  we  point  out  that  in  the  absence  of  a  finite  nonzero  density 
gradient,  certain  terms  vanish  from  the  cold-fluid  equations.  These  terms  are  also 
the  terms  which  are  responsible  for  the  resonant  diocotron  instability.  Thus  the 
box  profile  can  have  different  modes  than  the  ramped  profile  has. 


64 


D.  J.  Kaup  and  Gary  E.  Thomas 


For  now,  let  us  present  the  linear  equations  of  this  order.  They  are 


-iU'h  i  +  iknjilx  +  3,(nA.v)  =  0, 

(15a) 

A2 

~  '"A*  ~  ikh  ~  flA,  =  °’ 

(15b) 

-  iuJ5ly  -  +  Gu,,  =  0, 

(15c) 

[d2-k2)^x-  n1  =  0 

(15d) 

where  we  have  introduced 

kv0, 

(16) 

which  is  the  frequency  as  seen  by  the  streaming  electrons.  Since  the  anode  and 
cathode  are  conductors,  we  must  impose  the  boundary  conditions 

<^i(y  =  0)  -  0  =  4>1(y  =  /). 

(17) 

In  general,  this  overdetermines  the  solution  except  for  certain  values  (complex 
in  general)  of  u,  which  are  the  discrete  eigenfrequencies,  and  the  corresponding 
solution  is  called  the  eigenmode.  These  eigenvalues  depend  on  the  density  profile 
not  only  through  /i0  in  (15a),  but  also  through  v0  in  (16)  and  A2  in  (15b).  The 
standard  way  to  solve  for  the  eigenspectrum  has  been  to  reduce  (15)  to  a 
second-order  differential  equation  for  <£j,  as  in  Reference  [7],  However,  this 
introduces  a  false  singularity  wherever  w2  =  fi2.  This  can  be  avoided  by  reducing 
(15)  to  a  second-order  differential  system  for  the  velocities  instead.  Thus  let  us 
take 

=  iu}/k. 

(18a) 

%  =  p4>M, 

(18b) 

-  1  nA2  uu,  \ 

*»  )• 

(18c) 

“  iP'PQ, 

(18d) 

where  p  and  u  are  functions  of  y,  $  is  an  amplitude  and  is  only  a  function  of  the 
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slow  variables,  and 


If  we  also  define 


Q 

A 


B  = 


2A:«o  d,.n0 

QA2  dU'  ’ 

(19a) 

fl2  -  »2, 

(19b) 

*3,.»0  2k2n(]tf 
ftl42  ’ 

(19c) 

then  (15)  reduces  very  nicely  to  the  linear  second-order  system 

dyp  *=  Au,  (20a) 

dyu  =  Bp,  (20b) 

which  only  has  singularities  at  w,  =  0  and  ±  ft.  From  (17)  and  (18c),  one  can 
obtain  the  boundary  values  for  p  and  u. 


V.  Second-order  DC 

This  contains  the  quasilinear  theory  of  Davidson  (7)  and  also  has  been  covered  in 
[16].  The  equations  are  form-wise  identical  to  the  first-order  DC  equations,  (6), 
except  for  the  sources.  They  are 


«  Sg\  (21a) 

~  *  S£>,  (21b) 

Oog?  -  8$»  =  S2<°\  (21c) 

dfrf  =  Sg>,  (21d) 

where  as  before,  we  choose  n(30)  =  0  without  loss  of  generality.  Now,  the  sources 
are 

si0*’  =  -  ^r«0  -  *•*•)•  (22a) 

Sj?  “  ~  (5,V/„  +  c.c.),  (22b) 

Si,  "  ~  ~  t,o5xt;i°)  ~  "  (l'fcfiufii/  +  c-c-).  (22c) 

3?  -  -  a>0.  (22d) 


where  c.c.  stands  for  the  complex  conjugate  of  the  preceding  term. 
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Solution  is  now  straightforward.  Equations  (21  d)  and  (22d)  yield  a  unique 
solution  for  <f>(20>  upon  imposing  the  boundary  conditions 

*<°)(j,  =  0)  =  0  =  ** °>(y  =  l).  (23) 

Then  (21c)  and  (22c)  yield  the  solution  for  v (2J,  (21b)  and  (22b)  yield  u2”\  and 
(21a)  and  (22a),  upon  integration,  give  the  condition 

dTdy<t>0  ~  A2C2  =  -^(n?vly  +  c.c.)-^( +  c.c.) , 

where  C2  is  the  constant  of  integration.  Making  use  of  the  solution  (18),  this 
becomes 


where 


3rd  4>0  -  A2C2  =  Dd  n0, 


d  =  £z>* 


2y(V*)(p*p) 

*  •  (25b) 

and  y  is  the  imaginary  part  of  u.  (If  we  had  more  than  one  linear  mode  present, 
then  the  total  D  would  be  equal  to  a  sum  over  all  modes  [16]  of  the  individual 
£>*’ s.)  Differentiating  (24)  once  more  gives 

Srn0  +  C2dyn0  =  dy{Ddyn0),  (26) 


which  one  immediately  recognizes  as  a  diffusion  equation  for  the  density.  As 
before,  the  constant  C2  is  determined  by  integrating  (24)  from  the  anode  to  the 
cathode  under  fixed-voltage  conditions.  This  give 

/d [D{-dyn0)dy 

c‘  W7,  •  (27> 


which  for  dyn0  <  0  and  D  >  0  gives  C2  >  0. 

There  are  several  consequences  of  (26).  First,  although  a  general  arbitrary 
initial  condition  would  have  a  combination  of  modes  with  both  positive  and 
negative  values  of  y,  due  to  the  eJl"  term  in  (25)  one  expects  the  unstable  modes 
to  eventually  dominate,  thereby  driving  the  total  D  positive. 

Second,  as  a  function  of  y,  D  is  largest  wherever  the  magnitude  of  w,  or  A  is 
smallest.  The  values  of  y  where  this  occurs  are  the  positions  of  resonances.  We 
define  the  diocotron  resonance  to  be  where  the  real  part  of  ue  vanishes,  and  the 
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cyclotron  resonances  to  be  where  the  real  part  of  w,  is  equal  to  either  4  ft  or 
—  Q.  Thus  by  (25b),  diffusion  of  the  density  is  maximum  at  these  resonances. 
This  maximum  rate  of  diffusion  will  cause  plateaus  to  form  wherever  an  unstable 
mode  has  any  resonances  [16]. 

Third,  since  the  diffusion  rate  is  proportional  to  the  density  gradient,  a  box 
profile,  which  has  an  infinite  discontinuity,  is  unstable  against  diffusion,  and  will 
instantaneously  shift  toward  the  ramped  profile  shown  in  Figure  2(b). 

Fourth,  the  positive  value  of  C2  implies  that  there  should  be  a  net  upward  flux 
of  particles,  since  the  operator  dt  4  C2dy  will  tend  to  move  lines  of  constant  nQ 
upwards. 

Fifth,  this  diffusion  only  occurs  for  unstable  modes.  If  all  modes  are  stable, 
then  n0  does  not  diffuse  or  even  change,  and  the  background  remains  undis¬ 
turbed. 

Sixth,  any  slow  spatial  scales  in  ifr  (such  as  the  scale  x)  be  forced  into  the 
background  if  D  #  0.  This  follows  from  (24)  and  (25).  If  the  linear  amplitudes  \p 
are  functions  of  x»  then  so  is  D,  and  then  within  the  time  scale  of  t,  the  same 
slow  spatial  scales  will  appear  in  n0,  due  to  diffusion.  When  this  occurs,  the 
first-order  DC  equation,  (10),  becomes  excited,  and  the  density  then  starts 
changing  on  the  faster  time  scale  T. 

Seventh,  (24)  allows  equilibrium  profiles  to  exist.  As  mentioned  before,  in 
zeroth  and  first  order,  there  are  no  conditions  which  act  to  determine  the  profile 
shape  of  n0(y).  However,  here  in  second  order,  due  to  the  diffusion  mechanism, 
there  can  exist  profiles  where  n0  will  be  static.  To  see  this,  setting  dT<j>0  -  0  in  (24) 
gives 

=  A2(0)exp(c2/o"^y).  (28) 

Note  that  if  there  is  one  dominant  unstable  mode,  then  by  (27),  the  ratio  of  D/C2 
is  time-independent;  then  (28)  gives  a  unique  profile  shape,  which  presumably 
the  solution  of  (26)  would  relax  toward.  However,  determining  this  equilibrium 
profile  requires  knowing  D,  which  requires  the  eigenfunction  p(y)  [see  (25b)]. 
But  the  eigenfunctions  depend  on  the  background  profile,  so,  we  find  ourselves  in 
a  circular  situation.  How  one  can  break  out  of  this  is  not  known.  Although 
numerically  solving  the  coupled  system  of  equations  (20)  and  (24)  should 
converge  to  some  equilibrium  profile,  still  that  problem  remains  to  be  treated. 

Lastly,  I  wish  to  point  out  that  quasilinear  diffusion  is  not  the  only  important 
consequence  of  this  order.  Even  if  there  are  no  unstable  modes  (D  *  0)  and  n0  is 
stationary,  then  there  is  always  at  least  a  shift  in  the  shear  flow,  even  if  the 
background  has  no  slow  spatial  scales.  From  (18),  (20),  (21),  and  (22),  one  finds 


45 


n^0>-(^r+Mx) 


C,  + 


"o*xM 
fiA:  ) 


(29a) 


45 


-gdrdfa  +  +  c.c.)e2lr'. 


(29b) 
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[If  one  had  several  linear  modes  present,  then  the  last  terms  in  both  (29a)  and 
(29b)  would  be  replaced  by  the  appropriate  sums  over  k.]  These  last  terms  in  (29) 
are  always  present  regardless  whether  all  modes  are  stable  or  not. 

VI.  Second-order  fundamental 

The  main  consequence  of  this  order  is  to  determine  the  group  velocity.  Since  it 
may  be  useful  at  times  to  evaluate  these  higher  order  solutions  numerically,  we 
shall  briefly  discuss  several  methods  which  would  be  useful.  The  equations  of  this 
order  are  formwise  identical  to  the  first-order  fundamental  equations  (15),  except 
that  in  this  order  they  have  sources.  They  are 


where 


-iuen2  +  ikn#2x  +  9,(nefi2jf)  = 

(30a) 

-  iuj)2x  -  ik4>2  -  JfVly  = 

(30b) 

-  iwj} 2y  -  dy*2  +  U.v2x  =  S2">, 

(30c) 

dfo  2  -  k24>2  -  ft  2  =  S$, 

(30d) 

S*  =  -  9Tn t  -  dx(v0 nj  -  d^not^)  -  dy(htv^). 

(31a) 

5^  -  -  dT0lx  -  dx(voClx)+  dx4>t  - 

(31b) 

S$  =  -  drciy  -  v0dxvly  -  dy{v™vly). 

(31c) 

5<;>  -  -2ikdj>v 

(31  d) 

If  we  take  the  background  to  be  independent  of  the  slow  variables,  then  these 
sources  are  linear  in  dr\p  and  dy\p.  Thus  the  general  solution  will  be  of  the  form 


u2TdTyp  +  u2xdx4l 


V2*  “  k 

P2TdT'P  +  Pi  x'P 


(32a) 

(32b) 


Pit &  u2Tur 


P2X&  “2X“' 


dTi  +  v0dxi  ^uu>f/k  -  pA2/{ASl) 

*u — — —  +  («„*) - F - 

d  ip 

» 2  “  '(dT^PirQ  +  i{9x+)p2XQ  +  2pn 20-~ 


(32c) 
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where  p2T .  p2x.  uit-  an<^  u2x  sat'sfy  the  equations  (a=  27" or  2x) 


dypa  -  Aua  =  Ra, 

(33a) 

d,“a  ~  Vl>„  = 

(33b) 

+  u] 

R2T  =  2  ikn0p  ^2  , 

(34a) 

_  „  2  pnQu  uA 

R2x  “  voR2T  +  '  Ai 2  '  k  ’ 

(34b) 

.  .  /,  dyn0  4 k2n0tfur\ 

R2T~,p\kQU]A  AW  )• 

(34c) 

S  5  ■  >  ~2n0A7 

R2x  “  V0R2T  'pk  ^2q2 

(34d) 

Let  us  now  describe  how  the  solution  in  this  order  could  be  constructed.  The 
boundary  conditions  in  this  order  are  again 

•M-V*0)  “  0  «  $2(y  *=/).  (35) 

To  the  solutions  for  ( pa,ua )  we  can  always  add  an  arbitrary  amount  of  the 
homeogeneous  solution  (p,  u).  This  has  the  consequence  that  we  may  arbitrarily 
set  either  pa  or  u„  equal  to  zero  at  0.  At  y  =  0,  it  is  then  possible  to  set  the 
coefficients  of  both  dTt  and  dx^  in  $2,  (32c),  equal  to  zero,  by  judiciously 
choosing  the  initial  values  for  pa  and  u„.  One  then  integrates  (33)  forward  to  the 
anode.  The  solutions  for  pa  and  u„  are  now  unique.  In  general,  <^2  evaluated  at 
the  anode  will  be  of  the  form  frdT\p  +  /xdx^,  where  fT  and  fx  are  two 
constants.  Thus  (35)  then  demands 

dT\ p  4-  cdx^>  *  0,  (36) 

where  c-  /x//r  is  the  group  velocity.  The  major  consequence  of  this  order  is 
that  the  rf  amplitude  »|/  will  evolve  with  the  group  velocity  c  on  the  time  scale 
of  T **  et.  The  value  of  c  may  be  determined  by  the  above  procedure  of  inte¬ 
grating  (33). 

An  alternative  procedure  is  to  use  the  Wronskian  relationship.  From  (32a), 
(32b),  (33’  id  (20),  we  have 


dy(Auv2y  +  ikpvlM)  -  (dr*)(u/?2r-pK2r)  +  (dx*)(u/?2x-pfl2x)  (37) 
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Integrating  (37)  from  the  anode  to  the  cathode  results  in  (36)  and  the  initial  and 
final  values  of  v2x  and  v2y,  which  follow  from  the  boundary  conditions  (35). 

A  last  alternative  method  for  calculating  c  is  directly  from 


c 


du 
~dk  ' 


(38) 


If  one  differentiates  (15)  with  respect  to  k  and  x>  identifies 

h2  *=  —  idxdknl,  etc., 


and  sets 


(39) 


c(9x+)  -  -  dT+,  (40) 

then  one  obtains  exactly  (30)  and  (31).  Thus  the  solution  of  (30)  and  (31)  must 
give  a  value  for  c  which  is  equal  to  (38). 


VII.  Second-order  harmonic 

The  equations  in  this  order  are  again  formwise  identical  to  those  of  the  first-order 
harmonics,  except  that  w(/c)  is  replaced  by  2 u(k).  These  equations  are 


where 


1  +  2  ikn0vl£  +  dy(n  =  SQ, 

(41a) 

~  2'«,t -  2/A<&2'  -  ^  ujft  =  Sj2\ 

(41b) 

-  2 iajtffi  -  +  Slug  rn 

(41c) 

~  4 k7^22)  -  nf'  =  S<2\ 

(41d) 

“=  —  2/A«,t5Ijt  - 

(42a) 

SlV  =  -ik{Cu)2-  viyd)pu. 

(42b) 

Stf  «  -  ikvlx!)ly  -  6lyB/)iy, 

(42c) 

s$  «  o. 

(42d) 

The  only  significance  of  this  order  is  that  it  contributes  to  the  nonlinearity  in 
third  order,  and  to  fully  calculate  the  coefficient  of  nonlinearity,  it  is  necessary  to 
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obtain  the  solution  of  this  order.  The  general  solution  is  of  the  form 


of!  -  . 

(43a) 

eft  -  i’jf. 

(43b) 

,2/  P22&2  u22w»\  ,  , 2 P2B/A  — 

+>‘-l+\2Ajm  2*2)  *  2 *> 

V 

* 

(43c) 

(2)  .2/  Aknl  dy‘ 

"  ‘  ~  -  2 A 

'«)  +  *  !2»!".+2y“Mi^)  +  ?VJe;_ 

(43d) 

where 

n  .  1  a  Q  ,  a”o  at 

2utdyA  +  2kutA2dy\ 

5  k2\ 

A  A2) 

_  Us. 
^2 

- 

(44) 

A 2  =  S22  -4u2, 

(45a) 

4  k2- 

f) 

M^o)/(u»8)' 

-  &k2A2n0/(to2A2) 

(45b) 

a2 

^2 

1 

and  U22  and  p22  are  determined  by 

^ yPl2  ~  ^2U22  “  ^22* 

(46a) 

^>W22  ~  ^2  P22  = 

‘  &  22' 

(46b) 

where 

Ru  *=  4iufu2  -1 

-4/«#u/>d,(l//4)- 

~  ip^yiQ/A) 

-2iu,p 

.  (2  2o  /  ^ 

I2  >4  )' 

(47a) 

s  ,.,.2  u  +  P#yO/A) 

/?22  -  4/AA  u  Qa^ 

+  iktfp 

(2/4  2 

+  r(/t0- 

„  „  ,5,(*V/<J -*/><) 

-4--V 

(47b) 

(47b) 
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The  boundary  conditions  are  again 

=  0)  =  o  =  $2,(y  =  /).  (48) 

Integration  of  (46)  would  require  two  solutions  to  be  obtained.  The  boundary 
condition  at  the  cathode  will  determine  only  one  of  the  two  possible  initial  values 
for  P22^)  and  u22(0).  After  integration,  when  one  imposes  the  boundary  condi¬ 
tion  at  the  anode,  then  only  one  unique  linear  combination  of  these  two  solutions 
will  satisfy  both  boundary  conditions,  unless  an  eigenmode  also  exists  at  the 
harmonic  frequency  ( 2u )  and  wavevector  (2k).  In  general,  this  does  not  occur 
except  accidentally. 


VIII.  Third  order  DC 

In  the  case  that  there  are  no  unstable  linear  modes,  then  the  background  will 
remain  stationary  up  until  this  order,  which  contains  the  ponderomotive  force. 
The  equations  are 

3Anov3°y)  =  S3?,  (49  a) 

~  f  <  =  S'°x\  (49b) 

0»g?  -  dyW  =  S#\  (49c) 

W  =  S3?,  (49d) 

where 


•S3?  *  -  dr»o  -  3Anovu)  -  dx(utxnl+c.c) 

-  dy(n;6ly  +  c.c.+  n?02y+c.c.),  (50a) 

S3?  =  -  dr.v0  +  dTvg  +  -  6tflx) 


"  ~  (^VA  +  cc.)  -  (1 

»2y3,P  1*  +C-C.), 

(50b) 

S3?  -  -  aT<>  -  dTv?)  -  Mx<>  -  dy{ 

+  (ikv*j)lx  +  c.c.)  +  (tkvfyvlx  +  c.c 

•) 

(50c) 

S3?  -  -  d  v,0). 

(50d) 

The  Planar  Cold-Fluid  Magnetron  Equations  73 

As  before,  from  the  above,  one  obtains  the  secularly  condition  (upon  integra¬ 
tion) 


+  3xJ\e^xnl+c.c.)  dy  =  0,  (51) 


where  C3  is  the  constant  of  integration  and  is  to  be  determined  from  the  constant 
voltage  condition. 

Now  consider  (51)  when  there  are  no  unstable  linear  modes  and  the  back¬ 
ground  is  independent  of  the  slow  spatial  scale.  By  (29)  and  (18) 

d,dy* o  +  A JC3  «  Fr9r(M)  +  FMM),  (52) 


where 


Ft 


-  2wo  o 
“  QA  d 


2«o«’o 


2k/>Q 

A 


Referring  back  to  (36),  we  see  that  (52)  reduces  to 


dT.dy<p0  +  A2C3  -  (Fx-cFr)dx{ft), 


(53a) 


(53b) 


(54) 


which  gives,  upon  differentiating, 

dt.n0  +  C3*yn0-  9x(r*)-9,(Fx-cFTl  (55) 

The  constant  C3  is  determined  by  integrating  (52).  Thus 


C3 1'tfdy  «  [*,(*•*)]  jV*  "  cfr)  (56) 


At  the  moment,  it  is  not  possible  to  discuss  the  general  nature  of  the  functions 
Ft  and  Fx,  and  their  integrals,  as  in  (56).  However,  that  is  not  essential  for 


74 


D.  J.  Kaop  and  Gary  E.  Thomas 


describing  (he  nature  of  (55).  This  equation  describes  the  motion  of  the  back¬ 
ground  density  due  to  an  rf  wave.  It  occurs  on  the  time  scale  r'  (  =•  e3i)  and  is 
slower  than  all  other  effects  described  so  far.  What  is  important  here  is  the  fact 
that  the  evolution  of  n0  is  driven  by  dx(\p*\p).  This  is  the  ponderomotive  force 
(up  to  a  constant  factor).  In  one  region  it  will  drive  the  density  to  lower  values, 
while  in  another  region  where  £x(«£*i P)  has  the  opposite  sign,  it  will  force  the 
density  to  higher  values.  Of  course,  as  this  occurs  and  as  the  background  density 
undergoes  variations  with  respect  to  the  slow  spatial  scale,  the  first-order  DC 
equations  do  then  key  in  and  become  active.  Consequently,  it  is  of  little  use  to 
pursue  this  expansion  beyond  third  order,  because  at  this  order,  the  first-order 
DC  equations  will  always  become  keyed  in,  and  will  thereafter  dominate,  being 
of  a  lower  order. 

IX.  Third-order  fundamental 

This  order  will  give  us  the  nonlinear  Schrodinger  equation  (NLS)  [17-19].  The 
equations  are 


3  +  ‘kn0Plx  +  dy(nQ%)  = 

(57a) 

A2 

~  '“Ax  ~  ik<h  ~  ^  %■  =  Six’ 

(57b) 

-  iwji Jy  -  dfa  +  QZix  =  S<», 

(57c) 

(  3y  ~  ^  2  )  *#*3  —  ”3  =  S^ , 

(57d) 

wl.ere  the  sources  are 

Sii)  *  cdxn2  -  drnl  -  dx(v0n2)  -  dx{n0v2x) 

-  +  nfvg  +  v?xfi[2)) 

-  -Mvf; + + w; + ar,*?*). 

Six  ”  <^2*  -  dx(v0v2x)  -  ik ( v(2%x  +  v fj}{£) 

~  dJu  +  dj>2  -  vtydfa  -  o(2y  j}lx  ■ 

-  "  °lyd/>U’ 

Si?  =  cdxC2y  -  v03xS2y  -  drvly  -  2ikv*J)^y 

-ik{vWly-vgu;y)-vtxdxvW 

~  dy(v2y3ly  + 

-  -2 ikdj>2  -  Bfa. 
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Before  continuing,  we  shall  make  some  simplifying  assumptions.  First,  we 
assume  the  background  to  be  independent  of  the  slow  variables,  and  we  take  all 
modes  to  be  stable.  Hence  u  is  real.  Then  by  (9)  and  (29) 

<!  =  o  =  v'Z 

(59) 

We  shall  also  ignore  all  harmonic  terms  such  as  v(22*  etc.  We  further  simplify  (57) 
and  (58)  by  factoring  out  the  dispersive  part.  This  is  done  by  differentiating  the 
first-order  equations  (15)  twice  with  respect  to  both  k  and  x.  using  (38)-(40). 
Then,  upon  comparison  with  (58)  and  (59),  one  observes  that  the  third-order 
solution  is  simply 

h 3  =  -  \32dlhx  +  8n3,  etc., 

(60) 

where  the  6-quantities  satisfy 

-iU'Sn3  +  ikn0Sv3l  +  dy{n08u)y)  =  55^ , 

(61a) 

-  iU'SC 3,  -  ik  6<p}  -  jfSOj,.  =  6S3(]\ 

(61b) 

-  iuf8v3i  —  3v6<£ j  +  Sl8vix  =  SSjJ.’, 

(61c) 

(d2-k2)s^-srh  =  55';'. 

(61d) 

and  the  6-sources  are  now  only 

6S">  - 

(62a) 

WjJ*  =  “ 

(62b) 

8S£'  =  -9vu„ 

(62c) 

8S$  =  0, 

(62d) 

where  the  operator  2  is 

s-  »,  +  ifa>s;-£(0)a„’. 

(63) 

and  from  (29)  and  the  above 

(64) 
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All  the  remaining  nonlinear  terms  only  involve  which  is  the  second-order 
shift  in  the  DC  shear  flow. 

To  obtain  the  NLS,  we  first  obtain  the  equations  for  Sv3x  and  8uu..  From  (61) 
and  (62), 


a,(A663,)  +  ikA80ix  =  -  Q8C3a>,  (65a) 

dv(-ikS£ix)-ABSCh.  =  (65b) 

where 


=  a,  SSa»  +  ik  8S£'  +  -j-  5Sa*  +  (66a) 

8C<'>  =  (66b) 

Now,  as  a  consequence  of  (20)  and  (65),  we  have 
dt.(AuSCir  +  ikp &vJx) 

(67) 

Lastly,  integrating  (67)  from  the  anode  to  the  cathode  will  give  the  NLS.  One 
must  use  the  boundary  conditions  (17)  and 

fi^3(y;  =  0)  =  0  =  6<^j(>’  =  /)  (68) 

along  with  (61b).  One  then  reduces  (67)  to 

*va’  +  ?^d*)*  + ****’“  °’  (69) 

where  one  can  show  that 


(70a) 

(70b) 


r 


with  J/~  given  by 
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JY 


eL 

A2 


tody"  Q 
ku2 


4nl 

Au. 


2«0 

U 


2  »0(»t)P.»e  Y 
-J-(  —  +  T  ")• 


(70c) 


Now  we  have  N  and  N  given  entirely  in  terms  of  the  background  quantities  and 
the  first-order  solution. 

From  this  form,  one  could  numerically  evaluate  these  coefficients.  However, 
from  the  general  structure  it  is  possible  for  us  to  make  certain  genera!  statements. 
Consider  first  the  W-term  in  (70a).  There  we  see  that  the  integrand  is  propor¬ 
tional  to  dv> j„.  Thus  the  integrand  has  a  nonzero  contribution  only  along  a 
density  gradient.  Furthermore,  unless  u,  or  A  has  a  zero  inside  the  density 
gradient  (thereby  giving  a  logarithmic  contribution  to  the  integral),  the  quantity 
N  will  be  real  for  real  «. 

On  the  other  hand,  the  second  term,  A,  will  in  general  be  complex  if  A  or 
ever  vanishes  inside  the  plasma  sheath.  As  one  can  see  from  (70c),  there  are 
apparent  singular  terms  whenever  <o,  or  A  has  a  zero  anywhere  inside  the  sheath. 
Thus,  whenever  Jf  has  a  singular  point  inside  the  plasma  sheath,  a  complex 
value  for  N  will  in  general  result. 

When  this  is  the  case,  (69)  reduces  to  a  nonlinear  Schrodinger  equation  with  a 
complex  coefficient,  namely 


1  ci2v 


(71) 


where  T  (  =  N/A)  is  complex  in  general.  If  the  imaginary  part  of  T  were 
negative,  then  a  nonlinear  instability  would  occur  (while  if  the  imaginary  part 
were  positive,  a  nonlinear  stabilization  would  occur). 


X.  Summary 

As  we  have  detailed  here,  the  lowest-order  instability  is  the  first-order  DC.  which 
is  yet  to  be  analyzed.  However,  because  this  is  a  long-wavelength  instability,  it  is 
not  expected  to  be  of  any  importance  until  some  lower-order  instability  has 
keyed  it  in.  In  particular,  the  second-order  DC  may  key  it  in,  on  the  second-order 
lime  scale,  if  there  is  a  linear  instability.  And  even  if  nothing  else  is  unstable, 
then  at  least,  after  the  third-order  time  scale,  the  third-order  DC  equations 
(ponderomotive)  will  definitely  key  it  in.  However,  before  this  last  occurs,  one 
would  expect  the  third-order  fundamental  equations  to  go  unstable  by  the 
nonlinear  instability  mentioned  in  Section  IX.  This  instability  will  be  on  the 
second-order  time  scale,  and  thus  faster  than  that  of  the  third-order  DC. 
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It  is  shown  that  electrostatic  Vlasov-Poisson  perturbations  that  propagate 
parallel  to  the  magnetic  field  in  a  planar  magnetron  are  stable  for  both 
an  isotropic  and  also  for  a  particular  anisotropic  (Ty  =  3 Tx)  temperature 
distribution.  The  inhomogeneity  of  the  electron  density  is  fully  incorporated 
in  the  analysis.  The  proof  makes  use  of  only  the  dispersion  relation  of 
Trivelpiece-Gould  type,  without  actually  solving  the  eigenvalue  equation. 
These  results  suggest,  not  unexpectedly,  that  these  modes  should  be  stable  for 
all  such  anisotropic  velocity  distributions. 


1.  Introduction 

The  stability  of  electron  plasmas  in  crossed  fields  has  been  the  subject  of  wide 
investigation  in  connection  with  the  operation  of  such  devices  as  magnetrons, 
diodes  and  crossed-field  amplifiers  (Buneman,  Levy  &  Linson  1966;  Thomas 
1982;  Chernin  Si  Lau  1984 ;  Davidson  &  Chang  1984,  1985;  Ott  et  al.  1985).  The 
diocotron  instability  (Knauer  1966  ;  Davidson  1974),  among  others,  is  the  most 
familiar  instability  encountered  in  the  cross-field  configuration;  it  is  a 
longitudinal  (electrostatic)  instability  with  the  wave  vector  k  perpendicular  to 
the  direction  of  the  operating  magnetic  field  so  that  the  electron  drift  velocity 
v0  is  parallel  to  the  k  vector. 

Non-neutral  plasmas  in  the  crossed-field  configuration  are  inherently 
inhomogeneous  in  the  density  distribution,  and  the  velocity  distribution  is  often 
found  to  be  anisotropic,  so  that  the  plasmas  can  be  characterized  by  the 
different  temperatures  along  appropriate  directions  (Kaup,  Hansen  Si  Thomas 
1985).  A  non-Maxwellian  energy  distribution  of  the  plasma  particles  can  often 
give  rise  to  electrostatic  instabilities,  so  it  is  of  interest  to  study  the  possible 
instabilities  of  such  a  distribution.  The  purpose  of  this  work  is  to  investigate  the 
possibility  of  any  unstable  longitudinal  waves  propagating  parallel  to  the 
magnetic  field  in  an  inhomogeneous  and  anisotropic  electron  plasma  as  would 
exist  in  a  magnetron  or  .elated  microwave  device. 

Although  it  has  been  asserted  that  a  parallel  component  of  the  wave  vector 
tends  to  reduce  the  growth  rate,  thus  stabilizing  an  otherwise  unstable  wave 
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(Buneman  el  al.  1966;  Crawford  1967),  there  has  been  no  specific  analysis  of  the 
high-density  (o>*  &  flj),  inhomogeneous  and  anisotropic  regimes  that  one  would 
expect  to  exist  in  a  microwave  device.  Thus  there  is  a  lack  of  understanding  of 
how  a  parallel  component  would  tend  to  stabilize  an  unstable  wave,  particularly 
for  the  Vlasov-Poisson  system  and  for  typical  inhomogeneous  and  anisotropic 
distributions  that  could  occur  in  the  crossed-field  configuration. 

In  this  paper  wc  provide  such  an  analysis,  showing  that,  even  with  a 
moderately  strong  anisotropic  distribution,  instabilities  associated  with 
longitudinal  waves  propagating  parallel  to  the  magnetic  field  do  not  occur.  In 
particular,  we  find  that  certain  general  features  of  the  plasma  dispersion 
function  may  be  used  to  demonstrate  the  stability  of  these  modes. 

In  homogeneous  plasmas  the  propagation  of  waves  subject  to  spatial 
boundary  conditions  is  characterized  by  an  eigenvalue  equation.  Our  proof  only 
requires  the  use  of  an  effective  dispersion  relation  of  Trivelpiece-Gould  (1959) 
type,  without  actually  solving  the  eigenvalue  equation,  f 

This  paper  is  organized  as  follows.  In  §2  the  basic  equations  are  introduced 
and  the  equilibrium  distribution  function  is  discussed.  The  general  procedure 
leading  to  the  eigenvalue  equation  is  also  described.  In  §3  the  linearized 
solution  of  the  Vlasov  equation  is  obtained  for  an  anisotropic  equilibrium 
distribution  function  with  three  different  temperatures.  Differential  equations 
describing  the  propagation  of  the  electrostatic  waves  are  examined  in  various 
limiting  cases.  Finally,  we  prove  the  stability  of  the  wave,  subject  to  the 
boundary  conditions  of  a  smooth-bore  planar  magnetron,  for  the  isotropic  case 
and  then  also  for  a  moderately  strong  (Ty  =  37^)  anisotropic  case.  The  results 
for  these  two  different  cases  are  identical.  Considering  the  general  case,  we  then 
argue  that  these  results  strongly  suggest  that  this  mode  will  be  stable  for  all 
similar  velocity  distributions. 

2.  Basic  equations  and  particle  orbit  in  unperturbed  fields 

To  describe  our  electron  plasma  in  a  planar  magnetron  (see  figure  1),  we  use 
the  Vlasov  and  Poisson  equations 

|+v.V/— l(E+i,xB„).V./-0.  (1) 

V.E=— 4  rrejfdsv,  (2) 

where  the  notations  follow  standard  usage  in  plasma  physics.  The  magnetic 
field  is  taken  to  be  entirely  external  and  constant  in  time  and  space,  and  we  take 
the  perturbations  to  be  purely  electrostatic.  We  Bhall  designate  by  a  subscript 
zero  the  equilibrium  (zeroth-order)  quantities,  which  are  the  steady-state 
solutions  of  (1)  and  (2)  : 

(v- V)/o~~(Eo  +  “v  *  Bo)-^»/o  =  °,  (3) 

V.E0  =  -4t reJ/Qd*r.  (4) 

The  relative  directions  of  the  crossed  fields  E#  and  B0  are  illustrated  in  figure  1. 

t  Strictly  speaking,  we  only  prove  ‘spectral  stability’,  and  not  'linear  stability’.  For  a 
discussion  of  these  two  definitions  and  this  point,  see  Holm  el  al.  (1985). 
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Fiourb  1.  Geometry  and  electron  orbits  in  the  planar  magnetron. 


It  is  well  known  that  (3)  is  satisfied  by  any  f0  that  is  a  function  of  only  the 
constants  of  the  motion  of  the  particle  in  the  unperturbed  fields  E„  and  B0.  We 
have  the  following  three  constants  of  motion  (Kaup  ei  al.  1985) 


Vx  =  rn(vx-Qv), 

(5a) 

P,  =  mvz< 

(5b) 

w  *  {m(v*x  +  vl  +  vl)-e4 >0(y), 

(5c) 

where  ft  =  eBJmc  is  the  gyrofrequency  and  <J>0  is  the  electric  potential 
associated  with  the  equilibrium  electric  field  E0(E0  =  —  V4>0). 

Thus  /o(p,,p,,uj)  satisfies  (3),  and  the  equilibrium  density  n0  =  jf0d*v  is 
determined  by  (4)  once  the  functional  form  of f0(px,p„u>)  has  been  given. 

Now  (1)  and  (2)  are  linearized  with  respect  to  the  perturbed  quantities 
/'  =  /— /o  and  E'  —  E— E0  to  obtain 

(6) 

m 

V<D')  (7) 

Equation  (6)  can  be  solved  by  integration  along  the  unperturbed  orbit: 


J-+v.V-i.(E.+rxB.).V./'  =  - 
VM>'  =  4t7  ej/’d'v  (E'  =  — 


/(x,  v,  1)  =  -  ^  f  dt’  V'<t>(x',  t') .  VT./0(x',  v'), 

mJ-CD 


where  we  shall  henceforth  drop  the  primes  on  /'  and  <t>\  We  resort  to  normal¬ 
mode  analysis  for  (7),  seeking  the  solution  in  the  form 

<l>(x,<)  =  Q(y)eit‘~U4, 


which  is  appropriate  for  a  wave  propagating  in  the  direction  parallel  to  the 
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magnetic  field  in  a  medium  with  inhomogeneity  along  the  y-direction.  Then  (7) 
and  (8)  yield 

drO(y,)eUr£'_‘"T,  (9) 

where  we  have  substituted  the  general  functional  form  of  the  distribution 
function  /0(x',v')  **  fo{px,p,,u>)\  r  =  t'—t,  and  the  orbits  should  satisfy  the 
initial  conditions 

iv  =  y'-y  =  °> 

it  =  z'-z  =  0, 
v'  =  V 

In  general,  (9)  is  an  integral  equation  for  the  potential  <t>.  Instead  of  solving  an 
integral  equation,  we  choose  to  expand  4>(y')  about  the  position  y  : 

W  =  +  <“) 

This  is  certainly  a  good  approximation  if  the  orbit  size  (gyroradius  A)  is  much 
smaller  than  the  scale  length  of  the  inhomogeneity  (L  =  |V£’0/£,0|_1).  Kaup  el  al. 
(1985)  calculated  single-particle  orbits  in  a  spatially  varying  electric  field  E0(y) 
and  constant  magnetic  field  B#  by  a  singular  perturbation  method,  assuming 
e  =  AIL  to  be  small.  We  write  down  their  orbit,  neglecting  terms  0{e*): 

vi  —  v0+A(l  cos  (At +  $), 

v'y  =  — ^4Asin  (At  +  0), 

Ail 

it  =  v0T  +  —  lain  (Ar  +  0)  — sin0], 

iy  =  ,4[cos(At  +  0)  — cos$J], 

where  v0  =  eE0(y)/mil  is  the  drift  velocity  and  A  =  (fl*  —  Avein0/m)\  = 
[fi*— w*(y)]l  is  the  reduced  gyrofrequency  (Kaup  el  al.  1985;  Prasad,  Morales  & 
Fried  1985). 

In  (12)  the  velocity  components  at  r  =  0  are  solved  in  terms  of  the  constants 
of  integration  A  and  <f> : 

vx  =  v^  +  Acostf),  vv  =  —AA  sin^.  (13) 

Substituting  (11)  and  (12)  into  (9)  and  carrying  out  the  time  integration,  we 
obtain,  with  the  aid  of  (13), 
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where,  in  the  respective  lowest  orders, 


7°  IcV'-w ’ 


l  =  K-^o)1/  f _ 1  ,  1  \ 

1  n*  \lcvt—w  kvt—w±  A  kv,—(o±2A/ 

+^(fa^rfa^±2i)-  (15c) 

In  (15)  terms  with  double  signs  (±)  are  summed  over  the  two  signs. 

3.  Eigenvalue  equation  and  proof  of  stability 

To  proceed  further,  it  is  necessary  to  have  a  specific  equilibrium  distribution 
function  f9(px,pt,  w).  Kaup  el  al.  (1985)  devised  a  model  distribution  function  to 
model  an  electron  plasma  in  a  crossed-field  planar  magnetron : 

/o(P*.P*.“0  •'iVe-^’exp^-^^pljexp^-— (16) 

where  to,  pz  and  p,  are  the  constants  of  motion  defined  by  (5)  and  N  is  a 
constant.  In  terms  of  the  velocity  components  vx,  vy  and  vt,  (16)  reads 

Uvx>  Vw*)  =  ^(y)exp[-^rx-^£nyj  Jexp^-^rjjexp^-^pvj), 

(16') 

F(y )  -  exp  [fie<t>9{y)]  exp  ^  (fi -  y)  j .  (17) 

Then  the  unperturbed  density  is  calculated  as 

(181 

Substituting  (16)  into  (14)  and  carrying  out  the  necessary  velocity  integrals,  we 
obtain,  after  lengthy  algebra, 

-  {**  + 4*6*71,,  *[1+  &  Z(£o)]}4>  =  0,  (19) 

where 

P— 2£.Z<C.H-2[l±^]f„Z(£„). 


(20  a) 
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*-K.2(£.)-12[i±ij^5]£,.Z(£„). 

w±nA(m6\l 

^±n~  k  (2J  ”  °’1,2^ 


1  f®  e~(* 

Z(C,-sLm 

which  is  the  plasma  dispersion  function  (Fried  &  Conte  1961),  and 


(2u6) 

(20c) 

(21  o) 

(216) 


eE«  = _ !_l^o 

y  *  mil  m£lfi  n0  dy  ' 


(21c) 


which  is  the  Larmor  drift  velocity.  It  is  legitimate  to  neglect  mu*  in  (19)  as 
compared  with  t/y,  since  the  Larmor  drift  velocity  is  much  smaller  than 
the  thermal  velocity  (vT)  because  of  the  assumed  Bmallness  of  A/L(u!vT  = 
(vT/n0)dn0/dy  —  A/L  <  1).  This  neglect  is  even  better  justified  in  the  cold-fluid 
limit. 

We  shall  now  consider  a  few  limiting  cases  of  (19). 


(i)  Cold-plasma  limit  (w/k  >  vT) 

Expanding  the  plasma  dispersion  function  in  powers  of  (k*/w*)  v*T  and  taking 
the  limit  as  vr-*-0,  wc  can  reduce  (20)  to 


P  =  - 


2fiA* 


<J(w*-A1)’ 

2/?A*/ _ 1 _ 1  \ 

S  \w*  — A*  (o*— 4A*/’ 


R  = 


2/7A*  1 

6  w*  —  4A* 


Substitution  of  these  expressions  into  (19)  yields 

dw'p/dy  _ ^ ^ ^  ^ g,  =  0  (22) 

dy*  w*-Q  *+v*pdy  \  a,*)  v*-Cl*  (  ' 

In  obtaining  (22),  we  have  put  y  =  /?  and  A*  «  ft*.  Equation  (22)  could  have 
been  obtained  by  starting  with  the  cold  fluid  equations. 


(ii)  Massless  guiding-centre  limit 

If  the  electrons  are  treated  as  a  massless  guiding-centre  fluid  (m-*-0)  then  P,  Q 
and  R  all  vanish,  to  yield 

^-k*d>-int*n9S[l  +  ZoZ{Co)}*  =  0.  (23) 

This  equation  could  have  been  obtained  by  employing  the  drift  kinetic 
approximation,  in  which  the  electron  dynamics  in  the  plane  perpendicular  to 
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B0  are  suppressed.  We  can  also  see  that  the  drift  approximation  is  consistent 
with  the  approximation  of  very  high  frequency,  u>  >  A.  In  this  case  £±  as  and 
P,  Q  and  R  all  become  zero. 

(iii)  Low-frequency  wave:  A*  =  Q*— w*  >  at* 

When  A*  ►  w*  and  ft  &  £,  the  terms  indicated  by  _  in  the  expressions  for  P, 
Q  and  R  in  (20)  are  negligibly  small  compared  with  tc  £(£<>)•  Equation  (19)  then 
becomes 


dy*L  4  00  AVJ  fi*  dy  dy 

-{i*  +  4ffe*n0<ni  +  Co  £(£»)]}  *  =  0, 

which  can  be  cast  into  the  adjoint  form 

^(l-«u*r^j-{fc*  +  4JTe*n0[i  +  £0Z(Q]}(i-a<4r1<t>  =  0, 


where 


,  CqZ(Co) 
Q*v  ' 


4 

3/5/y+T 


Regardless  of  which  case  we  are  considering,  in  order  to  determine  <l>  uniquely 
from  one  of  the  above  differential  equations,  it  is  still  necessary  to  specify  two 
boundary  conditions.  We  shall  consider  the  two  parallel  planes  y  —  0  and 
y  —  l  to  be  perfect  conductors,  as  in  a  planar  magnetron.  Then  the  boundary 
conditions  are  that  the  Fourier  amplitude  <t>  must  vanish  at  the  two  boundary 
planes: 

C*(0)  *  6(/)  =  0.  (26) 

Next,  an  effective  dispersion  relation  of  Trivelpiece-Gould  type,  which  is  the 
solvability  condition  for  the  differential  equation  subject  to  the  boundary 
condition,  is  obtained  by  operating  with  fg<t>*dy  (*  denotes  the  complex 
conjugate)  on  the  differential  equation.  If  we  perform  this  operation  on  (24), 
first  assuming  v  m  1  (fi  =  y),  then  we  obtain 


where 


Z'  =  ^£l«_2ll+foZ(&))]. 

dQ o 


Putting  Z'  =  Z'r  +  iZ't,  where  Z’r  and  Z\  represent  respectively  the  real  and 
imaginary  parts  of  Z’,  (27)  separates  into  real  and  imaginary  parts: 


—  27re*»0  <J|4>| 


•)  =  o. 


Equation  (29)  indicates  that  the  solvability  condition  is  that  either  Z\  or  the 
quantity  in  parentheses  be  zero.  However,  the  latter  condition  makes  it 
impossible  to  satisfy  (28),  since  the  first  term  of  (28)  is  positive-definite  and  the 


last  term  would  have  to  vanish.  Therefore  the  solvability  condition  for  the 
system  (24)  subject  to  (26)  is 

£;  =  lmp^j  =  0.  (30) 

We  can  see  that  (30)  is  also  the  solvability  condition  for  the  system  (23)  and 
(26),  observing  that  (23)  is  formally  obtained  from  (24)  by  taking  the  limit  as 
Q-*-  oo.  If  an  exactly  parallel  procedure  is  carried  out  for  the  cold-fluid  equation 
(22)  then  the  solvability  condition  turns  out  to  be  Im  (o>)  =  0,  which  is  the  cold- 
fluid  limit  of  (30).  Thus  electrostatic  parallel  waves  in  the  cold-fluid 
approximation  are  stable  in  the  planar  magnetron. 

Returning  to  (30),  we  examine  the  possibility  that  it  could  be  satisfied  with 
Im  {(i))  >  0.  The  tabulation  of  the  plasma  dispersion  function  (Fried  &  Conte 
1961)  indicates  that  the  only  points  on  the  complex  ^  plane  that  satisfy  Z't  =  0 
with  Im  (w)  >  0  are  those  points  on  the  positive  Im  (£0)  axis.  But  the  value  of 
Z'r  on  the  positive  Im  (£0)  axis  is  restricted  by  —  2  <  Z'r  <  0.  However,  we  can 
immediately  see  that  with  —  2  <  Z'r  <  0,  (28)  cannot  be  satisfied,  since  it  can 
now  be  reduced  to  a  series  of  positive-definite  terms.  In  summary,  the 
eigenvalues  of  the  systems  (24)  (with  v  =  1)  and  (26)  Bhould  be  selected  from 
those  points  satisfying  Z'(  =  0  in  the  lower-half  £<,  plane.  These  eigenmodes  are 
therefore  only  Landau-damped. 

It  is  a  foregone  conclusion  that  the  eigenvalues  of  (23)  will  also  give  only 
damping  as  observed  in  the  solvability  condition  (28)  when  we  take  the  limit 
fi-*-  oo. 

Finally,  we  consider  one  case  of  anisotropic  temperature  by  assuming  v  =  2 
(y  =  3/?  or  Ty  —  $TZ)  in  (24).  The  analysis  of  the  case  for  a  general  non-integral 
value  of  v  is  much  more  difficult,  but  the  following  result  for  v  —  2  strongly 
suggests  that  the  above  does  extend  to  general  values  of  v.  Similar  algebra  to 
that  employed  above  yields 

a'1  ld<t>  * 

dy 

+  J  [4ne1n,5|2fii  — 4U1ojpar  — w*)  — =.  0,  (31) 


rf< i> 1 

J  [l-2o'rw*+w«(ar,-a*)]  —  dy 

+  J  {4nein0S[i—ar(otp  +  2a.rClt  —  2iltojtJ)(al  —  a1)] 
+  lc*{  1  -  ar  )}  |<t>|*  dy  =  0, 


where 


“r  2 n*(1+*Zr)’  a<  "  4Q*Z<‘ 


Again,  as  a  solvability  condition,  either  Z\  or  the  quantity  in  braces  in  (31) 
should  vanish.  If  we  choose  the  latter  condition  then  (31)  and  (32)  combine  to 
give 


«  d$>  *  .  i  /4n*  .  *\ 

—  +  (1:* -I- 4^e*n0 £)  <fy  +  (aj  +  af)  J  wi ~  ^  j<fy  =  0. 


(33) 
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In  our  present  theory  the  last  term  of  (33)  can  never  be  negative  since  the  ratio 
of  |<t»|*  to  \d<t>/dy\x  is  0(A  */L*).  Thus  (33)  cannot  be  satisfied,  and  one 

of  the  solvability  conditions  for  the  anisotropic  distribution  must  again  be 
Z\  =  0.  The  possibility  of  growing  eigenvalues  is  rejected  since  the  other  solv¬ 
ability  condition  (32)  with  a<  =  0  can  again  never  be  satisfied  if  —  2  <  Z'r  <  0, 
which  gives  ar  <  0.  In  conclusion,  even  for  a  distribution  function  with  aniso¬ 
tropic  temperatures  (»'  —  2),  the  eigenmodes  are  subject  to  only  Landau 
damping. 

In  the  general  case  of  non-integral  values  of  v,  it  is  certainly  true  that 

«  0  is  a  solution.  However,  it  cannot  readily  be  verified  that  this  iB  the  only 
possible  solution  for  the  imaginary  part  of  the  Trivelpiece-Gould  dispersion 
relation,  owing  to  the  occurrence  of  irrational  powers  in  (24).  Thus  we  cannot 
at  present  exclude  the  possibility  of  the  existence  of  some  other  unstable  mode 
where  a,  #=  0.  However,  we  do  not  expect  any  such  modes  to  occur,  since  no  such 
mode  was  present  for  the  v  =  2  case.  We  shall  assume  that  no  such  mode  exists. 
Thus  instability  could  only  occur  for  a(  =  0,  and  again  we  have  ar  <  0  and  the 
real  part  of  the  Trivelpiece-Gould  dispersion  relation  being  positive-definite. 
Thus  unstable  modes  with  at  —  0  cannot  exist,  and  the  eigenmode  must  be 
Landau-damped  whenever  at  =  0.  The  total  exclusion  of  the  possibility  of  any 
other  modes  ( a(  =#  0)  would  probably  require  an  analysis  with  numerical  or 
analytical  solutions  of  (24). 
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Summary.  We  study  the  evolution  of  weakly  nonlinear  hydro- 
dynamic  disturbances  on  a  static  cosmological  background,  pay¬ 
ing  special  attention  to  nonlinear  modulational  instabilities,  so- 
litons,  and  self-focusing.  Our  model  consists  of  two  distinct 
nonrelativistic  fluid  components  coupled  only  by  gravitation. 
The  two-dimensional  (cubic)  Nonlinear  Schroedinger  Equation 
(NLS)  is  found  to  govern  the  long-term  evolution  of  the  envelope 
of  weakly  nonlinear,  nearly  plane-symmetric,  almost  monochro¬ 
matic  acoustic  waves.  Nonlinear  modulational  instability  may 
even  arise  within  a  range  of  wavenumbers  for  which  both  modes 
of  the  linearized  theory  are  (Jeans-)  stable,  leading  to  the  pos¬ 
sibility  of  soliton  formation.  Extrapolated  to  a  realistic  expand¬ 
ing  universe,  this  result  suggests  that  nonlinear  modulational 
instability  might  “switch  on”  before  the  linear  (Jeans)  instability. 
Moreover  -  in  contrast  to  the  one-fluid  case  -  the  two-fluid 
system  studied  here  also  exhibits  a  violent  nonlinear  self-focusing 
instability  of  the  type  observed  in  experiments  with  optical 
beams.  Provided  certain  restrictions  on  the  wavenumber  and 
initial  conditions  are  satisfied,  self-focusing  leads  to  a  steep  rise  in 
the  density  contrast  at  certain  isolated  points  in  two  dimensions, 
corresponding  to  lines  in  three  dimensions.  (Of  course,  the  pre¬ 
sent  theory  can  only  follow  the  evolution  of  self-focusing 
singularities  until  the  condition  of  weak  nonlinearity  is  violated.) 
We  also  find  some  evidence  for  the  onset  of  nonlinear  saturation 
of  linear  Jeans  instabilities.  If  present,  saturation  would  imply 
that  nonlinear  instabilities  might  dominate,  at  least  under  certain 
circumstances.  In  that  case,  the  usual  picture  of  biased  galaxy 
formation  at  the  Jeans  mass  scale  would  have  to  be  regarded  as 
oversimplified.  Independently  of  our  particular  model,  the  pos¬ 
sible  existence  of  nonlinear  modulational  instability  implies  that 
caution  should  be  excercised  when  interpreting  the  results  of 
certain  “n-body”  numerical  simulations  used  by  various  re¬ 
searches  in  large-scale  structure  and  galaxy  formation:  Any 
numerical  method  which  (in  effect)  filters  out  the  high-wave- 
number  part  of  the  initial  fluctuation  spectrum  may  seriously 
underestimate  the  effects  of  resonant  phenomena  which  can 
transfer  power  from  high  to  low  wavenumbers. 
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1.  Introduction  and  model 

It  is  becoming  increasingly  difficult  to  reconcile  standard  models 
for  the  formation  of  large-scale  structures  and  galaxies  with 
upper  limits  on  the  small-scale  anisotropy  of  the  microwave 
background  (Uson  and  Wilkinson,  1985).  Recent  observations  of 
large-scale  structures  such  as  filaments  and  voids,  together  with 
the  task  of  explaining  fundamental  properties  of  galaxies  (such  as 
their  characteristic  masses),  have  led  recently  to  exotic  proposals, 
such  as  the  existence  of  cosmic  strings.  Without  taking  a  position 
on  the  merits  of  such  exotic  proposals,  we  investigate  in  this 
paper  the  possibility  that  one  can  do  without  them.  We  will  see 
that  acoustic  waves  in  a  self-gravitating  system  of  two  decoupled 
nonrelativistic  fluids  may  suffer  nonlinear  (modulational)  insta¬ 
bility  and  may  form  self-focusing  singularities,  depending  on  the 
relative  densities,  the  wavelength,  and  the  local  initial  conditions. 

The  predictions  made  by  any  linear  theory  of  fluctuations 
strongly  depend  on  assumptions  made  about  the  “initial"  fluc¬ 
tuations.  In  contrast,  nonlinear  systems  often  exhibit  qualitati¬ 
vely  similar  behaviour  for  large  classes  of  initial  conditions.  The 
fact  that  nonlinear  effects  near  recombination  are  presumably 
still  weak  does  not  necessarily  preclude  strong  qualitative  devia¬ 
tions  from  the  predictions  of  linear  theories,  provided  the  evolu¬ 
tion  time  is  sufficiently  long  for  the  nonlinear  effects  to  accu¬ 
mulate  (Tomita,  1967,  197);  Juszkiewicz.  1981;  Peebles,  1970: 
Liang,  1976,  1977;  Vishnaic,  1982).  By  means  of  a  multiple-time¬ 
scales  approach,  we  will  see  that  linear  dispersion  and  cubic 
nonlinearity  may  have  comparable  effects  over  sufficiently  long 
time  scales.  The  evolution  of  a  wave  envelope  with  "slow  "  trans¬ 
verse  dependence  will  be  seen  to  satisfy  the  well-known  Nonli¬ 
near  Schroedinger  equation  in  two  dimensions. 

A  nonlinear  mechanism  for  the  (occasional)  formation  of 
filamentary  structures  by  nonlinear  self-focusing  was  recently- 
suggested  by  one  of  us  (Kates,  1986).  It  was  erroneously  claimed 
that  self-focusing  occurs  (for  certain  range  of  the  wavenumber)  in 
a  model  consisting  of  a  single  fluid  on  a  static  cosmological 
background.  The  correct  coefficients  for  the  one-fluid  NLS  equa¬ 
tion  will  be  given  in  the  Appendix  of  this  paper.  As  it  turns  out.  if 
only  one  fluid  is  present,  modulational  instability  can  indeed 
occur,  as  claimed,  but  the  coefficients  of  the  resulting  NLS 
equation  never  admit  self-focusing  singularities,  because  the 
linear  dispersion  term  does  not  change  sign.  In  the  present  two- 
fluid  model,  the  linear  dispersion  term  can  take  both  signs. 

Our  interest  in  studying  a  system  containing  two  fluids  is  also 
motivated  by  strictly  astrophysical  considerations.  Recently,  con- 


10 


siderable  evidence  has  been  found  in  support  of  a  dark-matter 
hypothesis:  suppose  the  predominant  constituent  of  the  universe 
since  recombination  has  been  some  form  of  dark  matter  that  does 
not  interact  directly  with  ordinary  matter  and  in  particular  with 
electromagnetic  radiation,  except  through  the  gravitational 
potential.  If  one  also  supposes  that  “biasing"  occurs,  it  is  possible 
to  avoid  gross  contradictions  with  the  observations  (Frieman 
and  Will,  1982;  Davis,  1985;  Vilenkin,  1985;  Wilson  and  Silk, 
1981;  Wilson,  1983;  Efstathiou  and  Silk,  1983;  Bond  and 
Efstathiou,  1984;  Vittorio  and  Silk,  1984).  Moreover,  the  in¬ 
dependent  evidence  for  dark  matter  is  rather  compelling  (Rees, 
1984). 

The  most  plausible  candidates  for  the  dark  matter  can  be 
treated  as  collisionless.  If  at  a  some  epoch  both  the  dark  and  the 
ordinary  matter  are  nonrelativistic,  it  is  reasonable  to  describe 
the  mixture  by  a  Vlasov-Euler-Poisson  system.  Quadratic  non¬ 
linear  effects  in  such  a  model  were  considered  recently  by  one  of 
us  (Kates,  1987).  Disturbances  whose  characteristic  length  exce¬ 
eds  the  Jeans  length  by  a  sufficient  amount  were  found  to  satisfy 
to  good  accuracy  the  well-known  Kadomtsev-Petviashvili  equa¬ 
tion  of  type  I  (Kadomtsev  and  Petviashvili,  1970),  which  can  be 
solved  exactly  by  inverse  scattering  and  admits  two-dimensional 
lump  solutions.  These  lumps,  which  represent  regions  of  negative 
density  contrast,  suggest  themselves  as  possible  mechanisms  for 
cosmic  voids.  The  model  also  predicts  strong  biasing,  i.e.,  the 
density  contrast  of  the  ordinary  matter  in  lumps  is  much  larger 
than  that  of  the  dark  matter. 

Under  what  circumstances  might  a  model  consisting  of  two 
noninteracting  fluids  be  relevant  for  cosmology?  First,  suppose  at 
some  epoch  the  universe  can  be  thought  of  as  containing  pri¬ 
marily  dark  matter  and  one  kind  of  visible  matter,  for  example 
photons.  If  the  velocity  dispersion  of  the  dark  matter  is  suffici¬ 
ently  small,  then  it  is  reasonable  to  describe  it  as  a  pressureless 
fluid,  i.e.,  dust.  If  the  ordinary  matter  is  described  as  an  ideal  fluid 
with  pressure,  ther  the  mixture  is  an  important  special  case  of 
our  model. 

Second,  suppose  that  one  ignores  the  presence  of  dark  matter. 
Then  the  universe  may  be  thought  of  as  containing  two  noninter¬ 
acting  fluids:  the  baryons  and  the  photons  Because  of  the  dif¬ 
ferent  Z -dependences  (Z  =  redshift)  of  the  background  densities 
of  photons  and  baryons,  the  relative  proportion  evolves  with 
time.  As  we  shall  see  below,  the  relative  proportion  turns  out  to 
be  one  of  the  critical  parameters  for  determining  which  modes  are 
modulationally  stable,  which  are  modulationally  unstable,  and 
which  are  modulationally  unstable  with  self-focusing. 

We  begin  by  constructing  a  homogeneous,  isotropic  back¬ 
ground  (Newtonian)  cosmological  model  characterized  by  an 
expansion  parameter  a(r),  a  total  background  density  p0(r),  and  a 
cosmological  constant  A,  satisfying 

—  CPo(0  +  A/3  Jfl(f)  (1) 

j'lPolt)  n(f)J]»0  (2) 

(The  notation  and  terminology  follow  that  of  Peebles  (1980) 
unless  otherwise  stated.)  We  suppose  that  the  density  fluctuations 
are  small  compared  to  the  background  density  and  that  the 
peculiar  velocities  (that  is,  the  velocities  with  respect  to  the 
background  expansion)  are  small  compared  to  the  speed  of  light. 
If  the  material  is  strictly  nonrelativistic,  one  obtains  Newtonian 


d2a(t) 

-  as 

dt1  L 


fluctuation  equations  on  a  background  satisfying  an  expansion 
law  appropriate  for  nonrelativistic  material.  One  may  also  incor¬ 
porate  the  case  in  which  one  of  the  fluids  is  a  photon  fluid.  Of 
course,  in  that  case  the  expansion  law  (2)  for  the  background 
density  is  modified. 

We  will  concentrate  our  attention  on  two-dimensional  flows 
and  dependence  on  two  spatial  dimensions.  The  extension  to 
three  dimensions  should  be  evident. 

Imagining  that  the  material  consists  of  ordinary  and  dark 
matter  (of  course,  the  model  is  equally  applicable  to  any  two 
decoupled  fluid  components,  visible  or  not),  we  denote  the  exact 
density  of  the  ordinary  matter  by  n(x,  y,  r)p0(r)  and  the  exact 
density  of  dark  matter  by  A /(x,  y,  t)  •  p0(r).  (That  is,  n  and  ,Y  are 
the  relative  concentrations.)  The  background  concentrations  are 
denoted  by  n0(t)  and  N0(t),  respectively.  We  suppose  that  the 
pressure  gradients  are  strictly  proportional  to  the  density 
gradients. 

As  in  the  one-fluid  case  (Kates,  1986),  it  is  convenient  to 
measure  the  time  in  units  of  r0  =  (4jtGp0(r0))'  1J  and  the  spatial 
coordinates  in  units  of  a0,  the  expansion  factor  at  some  time  r0. 
(such  as  recombination).  Following  Peebles  (1980).  we  introduce 
a  modified  Newtonian  potential,  and  we  remove  the  dimensions 
with  a0  and  T0.  Similarly,  we  express  the  peculiar  velocities  in 
units  of  a0/T0  and  denote  these  scaled  peculiar  velocity  vectors 
by  [u(x,  y,  t),  c(x,  y,  t),  0]  and  [L/(x,  y,  r),  F(x,  y,  t),  0]  for  the 
ordinary  and  dark  matter,  respectively. 

In  this  paper,  we  proceed  as  if  a(f)  and  p0(r)  were  constants. 
These  conditions  can  be  achieved  self-consistently  by  choosing  a 
nonvanishing  cosmological  constant  A  so  that  Eqs.  ( 1  )-(2)  admit 
solutions  in  which  a(t)  and  p0(r)  are  independent  of  the  time 
(Einstein  static  universe).  In  the  Conclusions,  we  will  discuss  the 
degree  to  which  one  may  hope  that  important  features  of  the 
dynamics  presented  here  persist  when  one  allows  for  a  realistic 
expansion  law. 

The  system  to  be  studied  consists  of  mass-conservation  and 
Euler  equations  for  each  fluid  constituent,  together  with  a 
Poisson  equation  for  the  modified  potential,  whose  source  is  the 
sum  of  the  density  contrasts  of  the  ordinary  and  dark  matter.  (We 
assume  that  the  rotation  of  each  velocity  field  venishes.)  We  thus 
seek  functions  n(x,  y,  j),  N(x,  y,  r),  u(x,  y,  t),  v(x,  y.  r),  C(x,  y,  r). 


F(x,  y,  r),  and  <p(x,  y,  r)  satisfying 

f,('>)  +  cJ(nu)  +  cr(m)  =  0  (3) 

c[(u)  +  ur,u  +  i.c!>u  +  c:cx(log  n)  +  <(>,  =  0  (4) 

c'.(i')  +  ucxv  +  vd,v + cJ£,(log  n)  +  d>,  =  0  (5) 

c,(N)  +  e,(NU)  +  c,(NV)  =  0  (6) 

c,(V)+Ve,U+  VdtU  +  C!cx(log  b')  +  <f>l**0  (7) 

£>,( F)+  V£,V  +  Vd,V+  CJi?,(log  N )  +  <p,  =  0  (8 ) 

el<t>  +  d},<t>=(n-n0)  +  (N-N0)  (9) 


The  quantities  cJ  and  CJ  are  constants  representing  typical 
propagation  velocities  (squared)  in  these  units  and  are  propor¬ 
tional  to  the  temperatures  of  the  fluids  (idealized  as  idea)  gasses 
with  an  isothermal  equation  of  state;  the  final  results  would  differ 
only  slightly  for  other  equations  of  state).  The  two  fluids  are 
coupled  only  by  gravitation.  Note  that  the  background  corres¬ 
ponds  to  the  solution  N  =  N0,  n=*n0,  with  all  peculiar  velocity 
components  and  the  potential  vanishing.  In  a  realistic  back¬ 
ground,  of  course,  the  time  would  enter  Eqs.  (3)— (9)  explicitly. 
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X  Linear  theory 

Before  we  study  the  effects  of  nonlinearities,  it  is  useful  to  obtain 
the  dispersion  relation  describing  the  simple  linearized  theory  of 
disturbances  about  the  background: 

(S  +  A'0)  (R  +  n0)  =  n0N0  (10) 

Rsw2-c2k2  (11) 

Ssa)2-C2k2  (12) 

A  special  case  is  sketched  in  Fig.  1 . 

One  branch  (the  "upper"  branch)  of  the  hyperbola  defined  by 
(10)  always  passes  through  the  origin,  the  other  (“lower"  branch) 
passes  through  the  point  (  —  2n0,2N0).  The  allowed  portions  of 
;?ch  branch  are  of  course  those  lying  in  the  region  of  the  R-S 
plane  corresponding  to  k2  >0.  Restricting  the  discussion  to  these 
allowed  portions,  we  note  that  the  upper  branch  always  corres¬ 
ponds  to  modes  with  real  frequency.  The  lower  branch  contains 
(in  general)  a  stable  and  an  unstable  (imaginary-frequency)  part. 
[If  one  of  the  fluid  components  is  pressureless  (dust),  all  of  the 
lower  branch  represents  imaginary  frequencies  ]  The  linear  insta¬ 
bility  will  be  referred  to  below  as  the  Jeans  instability. 

Along  the  asymptotes  S=-A'0,  R=-n0  in  the  allowed 
regions,  the  dispersion  relation  approaches  that  of  a  single  fluid 
with  given  density  and  sound  speed. 

If  both  frequencies  corresponding  to  a  given  wavenumber  are 
real,  there  are  in  general  two  propagating  modes  with  distinct 
group  velocities.  We  will  be  especially  interested  in  those  instabi¬ 
lities  which  are  not  present  in  the  linear  theory  but  do  occur  when 
nonlinearities  are  taken  into  account.  Moreover,  as  we  will  see 
below,  it  is  possible  that  the  growth  of  linearly  (Jeans-)  unstable 
modes  may  saturate  or  slow  down  due  to  nonlinear  effects.  Thus, 
even  at  wavenumbers  where  the  linear  theory  predicts  Jeans 
instability,  the  dominant  instability  may  be  the  nonlinear  one. 

Of  particular  interest  is  the  high-mass-contrast  case:  Suppose 
that  one  of  the  fluids  has  a  very  low  concentration  compared  to 
the  other,  i.e.,  n0/N0<z  1  (see  Fig.  1).  If  c2>C 3,  the  asymptotes 


Fig.  1.  Sketch  of  dispersion  relation  (2. J )— (2.3).  (3.21)  for  the  case 
n0  «  0.1,  A'0  «  0.9,  C3/c3  =  0.7.  The  unshaded  area  designates  the  allowed 
region  where  k 1  >  0.  The  upper  branch  is  automatically  Jeans-stable 
(cu3  >  0).  The  lower  branch  contains  an  unstable  domain  ui2  <  0. 
Asymptotes  are  at  R  ■  —  n0.  5  «  -  S0.  The  region  shown  corresponds 
roughly  to  -  2  <  <  2.  —  2  <  S  <  2 


cross  in  the  unphysical  region  corresponding  to  k2  <  0  Each 
branch  remains  near  one  asymptote,  and  one  of  the  asymptotes  is 
close  to  the  origin.  The  lower  branch  satisfies  approximately 

co3=— 1  +  C2k2  (13) 

and  will  be  referred  to  below  (occasionally)  as  the  “principal" 
section,  because  it  resembles  the  one-fluid  dispersion  relation. 
The  upper  branch  satisfies  approximately 

w2  =  c2k2  (14) 

and  will  be  referred  to  below  as  the  "acoustic"  section,  because  it 
passes  through  the  origin  like  the  dispersion  relation  for  acoustic 
waves  in  the  absence  of  gravitation.  If  c2  <  C2,  the  asymptotes 
cross  in  the  physical  region  at  k2  =  C2  —  c2.  The  two  branches 
cannot  intersect,  of  course,  but  rather  interchange  roles.  In  this 
case,  the  term  “principal  section”  will  mean  either  branch  when  it 
is  near  the  S  =  —  A'0  asymptote,  while  the  term  "acoustic  section” 
will  mean  either  branch  when  it  is  near  the  R=  -n0  asymptote. 
Finally,  in  the  high-mass-contrast  case,  the  critical  Jeans  wave- 
number  kj  approaches  1/C  as  n0/A'0  approaches  zero. 

Thus,  the  linear  theory  predicts  that  the  solution  of 
Eqs.  (3>— (9)  for  generic  initial  conditions  will  be  dominated  at 
large  times  by  the  unstable  modes  of  longest  wavelength,  which 
in  the  high-mass-contrast  case  have  a  growth  rate  of  about  unity 
(in  units  of  the  characteristic  time  scale  T0).  The  basic  length  scale 
w’hich  is  singled  out  by  the  dynamics  (as  opposed  to  the  initial 
conditions)  is  the  Jeans  length. 

Since  the  Jeans  length  is  presumably  a  decreasing  function  of 
time  in  our  universe,  the  prediction  one  would  obtain  by  extra¬ 
polating  the  lineari:ed  version  of  our  model  would  seem  to  be 
that,  at  a  given  mass  scale,  growth  can  only  occur  after  a  well- 
defined  epoch.  As  we  saw  above,  the  largest  mass  scales  would 
become  unstable  first.  However,  there  is  no  justification  for 
expecting  the  linearized  theory  to  provide  an  accurate  appro¬ 
ximation  to  the  full  theory  over  arbitrarily  long  limes.  (More¬ 
over,  the  limitations  introduced  by  neglecting  the  expansion  and 
by  assuming  a  nonrelativistic  model  become  more  serious  at  long 
times  and  large  lengths).  In  what  follows,  we  will  investigate  the 
possible  role  of  nonlinearities.  As  we  shall  see.  at  a  given  mass  on 
length  scale,  it  is  possible  for  nonlinear  instability  to  switch  on 
"before"  linear  instability. 


In  a  weakly  nonlinear  system,  nonlinear  effects  can  become 
comparable  with  linear  ones  only  if  they  accumulate  over  long 
times  or  distances.  The  strength  of  the  nonlinearity  thus  intro¬ 
duces  new  length  and  time  scales  into  the  problem.  In  order  to 
discuss  weak,  nearly  plane-symmetric,  nearly  monochromatic 
fluctuations,  we  assume  that  the  physical  quantities  of  interest 
can  be  developed  in  asymptotic  expansions  depending  on  the 
coordinates  and  the  small  amplitude  parameter  £  in  a  prescribed 
way.  Writing  the  £-dependence  as 

(15) 

(16) 

r~£Jt'}  +  £3t  j+  . . .  (17) 

N~N0  +  tNi+ciN1  +  ciN>+  . .  .  (18) 


3.  Derivation  of  the  two-dimensional  nonlinear  Schroedinger 
equation  for  acoustic  disturbances 


n  -n0  +  cnl+c2n2  +  cini+  ... 
u  —  cu, +i2u2  + c*u2  +  ... 
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U  ~tUi+c2U2  +  c3Us  +  .  .  .  (19) 

K~£JFJ  +  £>F}+  ...  (20) 

<t>~t<t>]+t2<t>1  +  c1<t>i+  . . .,  (21) 

we  note  that  the  leading  deviation  of  all  quantities  from  their 
background  values  are  of  0(c)  or  higher. 

As  has  been  discussed  in  various  papers  and  text  books  (Van 
Dyke,  1964;  Cole,  1968)  on  two-timing  methods,  one  demands, 
not  only  that  the  expansions (15)-(2l) should  represent  pointwise 
valid  asymptotic  expansions,  but  also  that  the  validity  be  uni¬ 
form  over  some  possibly  r-dependent  domain  of  validity  which 
may  become  large  as  c  -*  0.  In  this  case,  the  domains  of  validity 
should  extend  over  distances  of  at  least  order  1  /c  (times  a 
wavelength)  and  over  time  intervals  of  at  least  order  \/c2  (times  a 
period)  at  some  fixed  position  with  respect  to  the  “wave  packet”. 
We  introduce  the  slowly  varying  spatial  coordinates 

Xscx  mi 


as  well  as  the  slow  and  slower  time  variables 


where  I  is  a  fixed  linear  operator  acting  on  the  Ok”}  parts  Qp  of 
the  expansions  (15)-(21)  and  S„  is  a  nonlinear  function  of  the 
terms  Q,,  Q} ,  Q,.,,  and  their  derivatives.  Resonances  occur 
whenever  Sf  contains  terms  which  are  not  orthogonal  to  the 
homogeneous  solutions  of  the  self-adjoint  operator  L  Reson¬ 
ances  would  lead  to  secular  growth  and  thus  nonuniformity  of 
the  asymptotic  expansions  (15)-(21).  In  order  to  obtain  uniformly 
valid  expansions,  we  employ  a  multiple-scales  technique,  which 
consists  in  choosing  the  functional  dependence  of  the  envelopes 
in  such  a  way  as  to  render  the  sources  S,  orthogonal  to  the 
homogeneous  solutions  of  L.  In  this  problem,  this  condition 
amounts  to  selecting  those  terms  which  are  proportional  to  a 
linearized  propagating  wave  solution  and  demanding  that  the 
sum  of  all  such  terms  vanish.  Note  however  that  the  nonresonant 
parts  of  Sp  need  not  vanish,  but  rather  will  drive  Qp. 

It  is  convenient  for  the  moment  to  carry  out  the  calculation 
with  the  y-dependence  suppressed.  As  expected,  at  0(c) ,  the 
solution  takes  the  form  of  the  linear  theory  in  one  dimension: 


(33) 


Tmct  (24) 

rsc2t.  (25) 

Derivatives  of  functions  of  the  slow  variables  are  evidently  down 
by  factors  of  £  or  c2.  Since  we  are  interested  in  the  evolution  of 
slowly  modulated  wave  trains,  the  dependence  of  all  quantities  of 
interest  on  the  fast  variables  (t,  x)  is  taken  into  account  by  rapidly 
varying  phase  functions,  whereas  the  dependence  on  slow  vari¬ 
ables  occurs  in  a  multiplicative  amplitude  or  envelope  function. 
(It  is  this  envelope  which  will  turn  out  to  satisfy  the  NLS 
equation  in  the  slow  variables.)  Note  that  the  y-dependence 
occurs  only  through  the  slow  variable  Y.  In  first  order,  the 
assumed  form  is  thus 


n1=m,e,'  +  m*e-i*  (26) 

+  (27) 

N,  =M1ei*  +  MTe_i'  (2g) 

1/,-v,  (29) 

4] -/i  (30) 


where  m,,  Mjt  pJ(  v,,/,  etc.  are  functions  of  (X,  Y,  T,  t).  It  is 
convenient  to  introduce  the  notation 

,  d0  ee 

km—  um—  (31) 

cx  dt  '  ' 

for  the  derivatives  of  the  rapidly  varying  phase.  Note  that  to  is  the 
negative  of  the  physical  frequency.  [In  principle,  the  frequency 
could  be  large  compared  to  unity,  as  long  as  it  is  restricted  to 
satisfy  a>24  1/e,  in  order  that  the  higher-order  terms  be  truly 
decreasing  in  order  of  magnitude  (Kates,  1986).  However,  no 
generality  is  lost  by  treating  the  wavenumber  and  frequency  as  if 
they  were  just  of  order  unity,  as  will  be  done  here  ] 

Our  procedure  is  to  collect  terms  order  by  order,  paying 
special  attention  to  resonances:  Now,  at  each  order  cr,  the  system 
assumes  the  form 

LQ,mS,>  (32) 


c o  m. 


k  n0 

(34) 

w  /, 

2  /-'I’ 

or  —  C‘k‘ 

(35) 

w  M , 

'*  ~k~K^ 

(36) 

,  _  n0  Nc 

ta2—c:k2  w2  —  C2k2 

(37) 

where  the  dependence  of  /,  on  the  slow  variables  is  as  yet 
undetermined.  (Note  that  (37)  is  equivalent  to  (10)-(12|.) 

At  0(c2 ).  the  solution  contains  contributions  at  the  funda¬ 
mental  and  iwice  the  fundamental: 

*=  (m'j1  V  +  c.c.) +  (m'/'e2'*  +  c.c.) 

(381 

N  a  =  'e*  +  c.c.)  +  (Mf k2‘e  +  c.c.) 

(39) 

ttj  =  (Pi  V*  +  c.c.)  +  (ff'e2'*  +  c.c.) 

(40) 

U:  -  (v-y  V*  +  c.c.)  +  (v';V*  +  c.c.) 

(41) 

0 J  =  (/:  v**  +  C.c.)  +  ( /'/ V ***  +  c.c.). 

(42) 

Without  loss  of  generality,  we  set 

/a*-0. 

(43) 

The  solutions  for  the  terms  proportional  to  e2",  e~2il 
form 

’  take  the 

l*j'm  “^['"iJ'/no-K/n0)1] 

(44) 

v'/'« -^[M'/'/No-tM, /No)*] 

(45) 

*#'_  *:/,2ji  ir,  2w2  ir«,T 
fl0  w2 ~c2k2  +  2l_  1  w2-c2k2  Jl_rt0  J 

(46) 

m'/>  k2r>'  ir  2w2  -jfM.i3 

No  (uJ-CJkJ+2L  <u3-CjH!JLn;J 

(47) 
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■-Mu 

6k3L  w3-c3k3JUaJ 

_ IF— T 

6k3|_  V-C3k3JLwJ  ' 


The  terms  proportional  to  e ",  e'M  yield  several  algebraic  rela¬ 
tions  as  well  as  the  first  secular  condition: 

4 n-  ~T—  +  -W[-a>c,  +  kcr]m,[X,  T,  r.i]  (49) 


(w2-c2k2): 


;[-wdx+kcTy,[X,  Y.T.  t] 


•i1'-—' +^[-«cjr+kcr]M,[x,  y,  r,  i]  (51) 

Mi" = [  -  cue, + *cT]/,  [X,  y,  r.  o  (52) 

^(/,)+— t«f*,-«r]/.-a  (53* 

w 

where  we  define  for  convenience 
_  n0co2L~ 2  N0 w2L~2 

L  (a)2—  c2k2)L  (to3  —  C2k2)L  ^  ^ 

The  secular  condition  (53)  is  a  property  of  the  linear  theory  and 
implies  that,  to  a  first  approximation,  the  amplitude  information 
propagates  along  the  integral  curves  of  the  vector  field 

(55> 

-(1 +Dj)— ,  (56) 

k 

where  vl  =  d(  —  w)/dk  is  the  group  velocity  of  the  wave  packet. 

At  0(c3),  the  equations  split  into  part'  proportional  to 
e12'*,  and  e *’*.  The  parts  of  the  0{c3)  solutions  proportional  to 
et3it  and  ex2,t  exert  no  influence  on  the  secular  terms  at  this 
order  and  will  play  no  role  in  what  follows.  Let  us  denote  the 
parts  proportional  to  e1"  by  m'j",  M'3n,  p'j ’,  v'j", /V1,  etc.  as 
before.  We  collect  only  those  terms  in  the  equations  proportional 
to  e±u.  Using  the  mass  conservation,  Euler,  and  Poisson 
Eqs.  (3)— (9)  and  setting  (without  loss  of  generality)  /"'sO,  we 
obtain  the  relations 

»o^V  ’  “  “  ' +  i  +  ^ r«j  * + '] 

(57) 

J  °U  2  n30  JV**3 

to3-c3k3|_  L"o  "o  J 

+ft>  Tc*brJ 

- ’+ ^|]  +  ia,[^[~]  +  ^[' +  ’  J 


fV!>"  =  -  7  MV 1 + kc,  M ,  +  cT  wy 1  +  A'orx 'V  ’] 

K  K 

(59) 

2  A’o  J  * 

+^_r,:v.v.,_4^_vy1+^vTii 

w3-C3*3[_  L^o  ,w0  JJ 

jT M2M*  io>2  rwyn 

“  L  w§  “  J  k  e*l  No  J 

-  i/c[c7.vyl  +  ctv1]  +  +  +  cj.v'j1 1 

(60) 

-^Sr/.+mV'+My-O.  (61) 

Substituting  (37)  and  (59)  into  (61)  and  using  (33),  (34),  (52)  and 
(54),  we  obtain  the  following  expressions,  without  the  term  in¬ 
volving  Cydrifi )  in  braces. 

2 ik2  k6  T  12  1 

- —  2Dj  +  £>4  — -Dj  — -DjD3  — -[Dj]3 

w  u  L  2  3  6 

-^[D»]Jl/T(/i)J  +  -T[i)J-4®»3 

3  J  or 

,  2  D} 

X  [  -  wcx  +  kcTyf,  +  —  [wr*  -  kcT’] r,/,  -  <*(-,/,  + 

Cl) 

{-[1  +r»2]crcr(/,)}  =  0.  (62) 

Using  the  secular  condition  (53)  we  find  (again  without  the  term 
in  braces) 

W  CT 

*6T  12  1  ,  2  ,1 

+  42  i  +  £>4“2Dl_3D2D3_6[£)j]  3tDj] 
r-3(DJ)3  +  Dj-4Dj'|f3/,  f  c3/,] 

T  ~W  (63) 

Without  the  term  in  braces,  Eq.  (63)  is  the  one-dimensional 
nonlinear  Schroedinger  equation.  However,  it  is  possible  to 
generalize  the  above  calculation  in  order  to  take  into  account  the 
y-dependencc  as  follows:  The  linear  terms  in  (63)  are  those  which 
would  be  obtained  by  applying  the  present  singular  perturbation 
method  to  the  linearized  version  of  (3)— (9).  It  can  be  easily 
verified  that  the  nonlinear  term  is  unaffected  by  the  inclusion  of 
y-derivatives  and  y-components  of  velocity.  We  have  therefore 
carried  out  a  second  calculation  of  the  linear  coefficients,  includ¬ 
ing  y-derivatives,  by  constructing  the  linear  operator  associated 
with  the  dispersion  relation  (10)  and  expanding  in  terms  of  slow 
derivatives  using  the  computer  algebra  system  MUMATH.  This 
procedure  confirms  the  coefficients  obtained  above  and  provides 
the  y-coefficients  in  braces  in  (62)- (63). 

Finally,  letting 

k2 

F-—f  (64) 

ti r 

{-kX  (65) 
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q  *  k  F  (66) 

f«on,  (67) 

we  obtain  the  two-dimensional  nonlinear  Schroedinger  equation 
in  the  form 


icF  c2F  c2F 

— +2WlFmF2]  +  D—  +  £—  =  0,  (68) 

of  ci 1  Ct) 1 

where 

^Sii-[2D5  +  D4-iD3-^JD3-I[D2]^^[D3]>]  (69) 

D*^[D>-HD>V-*D>1  (70) 

£M-j[l  +  Da]/D2.  (71) 


The  result  of  our  multiple  time-scales  approach  is  thus  the 
following:  Provided  F  and  its  complex  conjugate  F*  are  chosen 
in  such  a  way  as  to  satisfy  Eq.  (68)  and  its  conjugate,  no  secular 
terms  will  arise  in  the  solutions  up  to  and  including  0(e3),  as  is 
necessary  if  our  expansions  (1 5)— (2 1 )  are  to  represent  uniformly 
valid  asymptotic  approximations  (over  the  distance  and  time 
scales  considered)  to  some  exact  solutions  of  the  original  problem 
(3)-(9)  corresponding  to  particular  initial  value  data.  [Eq.  (68)  is 
a  necessary  condition;  whether  it  is  sufficient  for  uniformity 
remains  an  open  question  ] 


4.  Qualitative  properties  of  the  two-dimensional  NLS  equation 


The  degree  to  which  the  nonlinear  theory  differs  qualitatively 
from  the  linear  one  depends  on  the  relative  signs  of  the  co¬ 
efficients  D,  £,  and  W'.  Since  any  solution  of  the  one-dimensional 
NLS  equation  also  solves  the  two-dimensional  version,  it  is 
convenient  to  review  a  few  of  the  properties  of  the  one-dimen¬ 
sional  NLS  equation  expressed  in  the  form 


icF  ,  c2F 

— +  2K'F*F3  +  D— =  0. 

fT  Cf 1 


(72) 


If  W and  D  take  different  signs,  then  solitons  do  not  form  and  the 
underlying  acoustic  disturbance  is  called  “modulationally  stable” 
(i.e.,  small  modulations  of  a  plane  wave  are  always  stable),  and 
the  predictions  of  the  nonlinear  theory  are  not  much  different 
from  those  of  the  linearized  theory.  If  N  and  D  take  the  same  sign, 
then  the  acoustic  disturbance  is  said  to  be  "modulationally 
unstable”  (i.e.,  there  exist  unstable  perturbation  modes  of  a  plane 
wave),  and  soliton  formation  is  possible. 

Now,  the  one-dimensional  NLS  should  give  a  good  descrip¬ 
tion  of  the  dynamics  during  some  time  period  if  either  one  of  the 
terms  of  (72)  involving  spatial  derivatives  is  small  compared  to 
the  other.  Thus,  if  the  transverse  dependence  is  slow,  even  with 
respect  to  F,  then  solitons  can  be  expected  if  W  D>0.  On  the 
other  hand,  it  is  also  possible  in  principle  for  the  disturbances  to 
be  so  coherent  that  the  X -variations  are  much  smaller  than  the 
F-ones,  and  in  this  case  "transverse”  solitons  would  be  possible  if 
♦FcO.  (Note  that  £< 0  always  holds  in  this  system.) 

Although  of  course  solutions  of  the  one-dimensional  NLS 
equation  also  satisfy  the  two-dimensional  version,  they  are  un¬ 
stable  over  very  long  time  scales  with  respect  to  slowly  varying  or 
small  perturbations  in  the  second  spatial  direction  (Ablowitz  and 


Segur,  1981).  Therefore,  solitons  should  not  be  expected  to  persist 
indefinitely,  even  if  they  do  form.  Unfortunately,  the  general 
solution  of  (68)  is  not  yet  known.  However,  solutions  of  (68) 
satisfy  a  sequence  of  integral  identities,  the  first  few  of  which  are 


■■\p 

II 

© 

(73) 

CT 

f£i=o. 

(74) 

cf 

c2J2 

-^r=8J„  - 

ct1 

(75) 

where 

^([f-f^ 

(76) 

J,  =  JJ [£»  ■  [F{]2  +  £  ■  [F,]3  -  WF2F*2yidrt 

(77) 

Jj«|J[D{J  +  £q*]F*Frf{</ij- 

(78) 

If  the  coefficients  D,  E  and  W  are  all  of  the  same  sign  (say 
positive)  then,  for  certain  nonsingular  initial  data.  Eq.  (68)  exhi¬ 
bits  self-focusing  singularities  in  a  finite  time  (Zakharov  and 
Synakh,  1976),  as  one  can  make  plausible  by  the  following 
argument:  Suppose  the  solution  for  F  were  regular  at  all  times. 
For  the  case  D  £>0,  the  integrand  of  J}  is  positive  definite  On 
the  other  hand,  if  the  initial  data  has  the  property  that  the 
integral  is  negative,  then  by  virtue  of  (74)  it  will  remain  so 
However,  Eq.  (75)  implies  that  a  time  f0  exists  such  that  J 2<0. 
leading  to  a  contradiction.  The  nature  of  self-focusing  singu¬ 
larities  has  been  investigated  theoretically  (Talanov,  1965; 
Zakharov  and  Synakh.  1976;  Newell,  1978).  and  their  occurrence 
has  been  observed  in  nonlinear  optics  (Chiao.  1964,  Kelley.  1965). 

Of  course,  only  the  onset  of  self-focusing  is  a  true  prediction  of 
the  theory  presented  here:  If  the  density  contrasts  become  large 
compared  to  unity,  then  the  problem  enters  the  fully  nonlinear 
regime  and  the  expansions  assumed  here  are  no  longer  appro¬ 
priate.  The  essential  point,  however,  is  that  the  predictions  of  the 
nonlinear  theory  differ  qualitatively  from  those  of  the  linear 
theory:  Even  at  wavenumbers  where  the  linear  Jeans  instability  is 
not  present,  the  nonlinear  theory  predicts  the  possibility  of  soliton 
production  or  the  possible  occurrence  of  occasional  dramatic 
pointlike  increases  in  the  density. 


5.  Regimes  of  stability,  modulational  instability,  and  self-focusing 

It  is  routine  to  work  out  the  symptotic  values  of  the  coefficients 
W7,  D,  and  £  along  both  branches  of  the  dispersion  relation  in  the 
general  case: 

Along  the  principal  section  (lower  branch),  asymptotically  in 


the  limit  a>3-*oc. 

To)3  T"’ 

J  +0[(«3/n0)L-3] 

(79) 

£  =  —  1  /2  +  0[(cu3/n0)]  " 3 

(80) 

+0[(wJ/n°rJ] 

(81) 
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w°~\  Or  J+0[(a,2/"o)2}  (82) 

The  asymptotic  values  in  the  limit  a>J-»oc  along  the  acoustic 
section  (upper  branch),  can  be  obtained  from  (79)— (82)  by  replac¬ 
ing  n0  with  N0.  Thus,  at  sufficiently  high  frequencies  along  both 
branches,  D  is  always  positive,  W  is  always  negative,  there  is  no 
self-focusing,  and  in  fact  acoustic  disturbances  are  modulation- 
ally  stable  with  respect  to  the  longitudinal  direction.  Since  £  is 
negative,  the  only  possibility  for  soliton  formation  predicted  here 
is  the  case  where  the  acoustic  disturbance  is  very  coherent.  Then, 
as  described  above,  the  one-dimensional  NLS  equation  in  the  Y 
and  r  variables  admits  soliton  solutions. 

In  the  limit  tuJ-»0  along  the  upper  branch  (acoustic  section). 


the  values  approach 

DL  =  dLk-2  +  0(  1) 

(83) 

D=-^-k2  +  0(k*) 

2«: 

(84) 

£=  --  +  0(AJ) 

2 

(85) 

W=-±[2di  +  d2¥d-2'k-2, 

24 

(86) 

where 

j  (n0C2  —  K 0c2)L~ 1  r  .  k,  ,fl  (-1)1] 

dL=  ( cwy  K+Ao]U  M  J 

(87) 

Thus,  all  three  coefficients  are  negative  and  self-focusing  can 
occur.  Note  that  the  coefficient  W  diverges  like  k~2. 

Finally,  as  aj2—0  along  the  lower  branch  (principal  section), 
we  obtain 


Dl  =  u>2l~2(—  \)LdL  +  0(w2L) 

(88) 

0=-=^— +  0(cu-J) 

2  (d2)2w* 

(89) 

1 

£  - = — t-  0(1 ) 

2w2(d}) 

(90) 

JF=^l[3d,-(JJ)J]  +  O(cu4)>0, 

2Aa2 

(91) 

where 

Wq  N0 

^+c51 

^Lm(  1)\-  w  1L 

r>0  N0 

(92) 

.c1  +  C2_ 

Self-focusing  does  not  occur  here,  but  since  W  D>0,  the  system 
is  modulationally  unstable  with  respect  to  the  longitudinal  direc¬ 
tion,  and  solitons  may  form. 

For  reasons  described  in  Sect.  1,  it  is  also  of  interest  to  study 
the  high-mass-contrast  case  (n0/N0  <  1 )  in  detail.  In  that  case,  we 
can  develop  the  coefficients  in  powers  of  n0/N0  along  either 
section.  The  cross-over  region  (which  occurs  only  if  c 2  <  CJ) 
requires  special  care. 

Along  the  principal  section  (which  by  definition  means  out¬ 
side  the  cross-over  region),  the  frequency  is  approximately 


(o2  =  -  1  +  C2k2  +  n0A  +  0|>o/No):]. 

(93) 

where 

1 

(94) 

The  coefficients  are  then  given  approximately  by 

DL=(  —  \)Lw2L~ 2  +  ... 

(95) 

lF=[wJ(l  +  5«J-8w4)/12]+  ... 

(96) 

D  =  c2k2/{ 4a>4)+  . .  . 

(97) 

£=  -c2k2/(  2w2)+  .... 

(98) 

Along  the  acoustic  section,  the  frequency  satisfies  approximately 

w3  =  c2k2  +  n06+  .  .  . 

(99) 

(C2-c2)k2 

1-(C2-c3)k2' 

(100) 

and  the  coefficients  are  approximately 

DL  =  (w2/n0)L~’t 5L+  .  . . 

(101) 

N  =  [wJ/(n0i5)]J[2, 3—  l/(6kJ(CJ  —  c1))]  +  .  .  . 

(102) 

D=  —  [(n0<5J/(2w2)][l  -t-3/(i::(C2  — c2))]+  .  .  . 

(103) 

1 

£=— +  .... 

2 

(104) 

The  limiting  values  along  both  branches  for  large  and  small 
frequency  agree  with  those  found  above  in  the  general  case. 

In  the  cross-over  region,  both  hyperbolae  are  near  the  asy¬ 
mptotes,  and  a  different  expansion  is  called  for.  To  leading 


‘•order”  in  (n0/N0), 

wJ  =  cJAJ  +  (n0)'  Jr+  .  .  .  (105) 

(C2  —  c2)k2=  1  +2(«0)‘  2d  +  .  .  .  (106) 

r  =  d±(l+d:)‘  5,  (107) 

with  the  upper  (lower)  sign  corresponding  to  the  upper  (lower) 
branch,  (see  Fig.  1).  The  coefficients  become 

Dj  =  w:(l/r:  +  l)-l-  .  .  .  (108) 

Oi  =  (wJ/n0)1’” '(tto)1  2/TL+  ...(£>  3)  (109) 

N  =  w6/(2r3(rJ+l)nS'J)+  .  .  .  (110) 

D=  —  2r,/(cuJ(n0),/J(rJ+l),)+  ....  (Ill) 


In  the  one-fluid  limit,  which  is  discussed  in  the  Appendix,  the 
three  coefficients  W,  D,  and  £  are  never  all  of  the  same  sign  and 
self-focusing  never  occurs. 

It  is  also  of  interest  to  examine  the  effects  of  nonlinearities  on 
the  linearly  (Jeans-)  unstable  modes  where  a>J<0.  Now  the 
interpretation  of  a  “secular”  condition  changes  of  course  when 
the  underlying  disturbance  has  imaginary  frequency.  However,  it 
is  still  necessary  to  enforce  the  secular  conditions  (55)  and  (63)  in 
order  to  keep  the  higher-order  terms  in  ( 1 5)— (2 1 )  from  growing 
larger  than  the  lower-order  ones  (e  g.,  the  function  c3te"  would 
eventually  grow  larger  than  te'").  Taking  wsiy,  and  noting  that 

y*-l -c2k2,  (112) 

we  examine  the  implications  of  the  conditions  (55)  and  (63)  as 


k-> 0,  which  is  where  the  linear  instability  grows  the  fastest.  The  The  typical  time  for  the  nonlinearity  to  have  an  effect  depends 
conditions  become  on  W  and  is  roughly 


if  t/i  “  0 


(113) 


r  crf,  +  2kV’*/T  f\  +  -  C3  [r }  +  ciU  =  0. 


(114) 


Equation  (1 14)  is  a  (nonlinear)  diffusion  equation,  and  gradients 
tend  to  get  smoothed  out  with  time.  The  nonlinear  term  would 
then  have  a  tendency  to  cause  /,  to  saturate.  Ignoring  spatial 
derivatives,  we  find  that  the  quantities  /f/i  and  the  potential  4> 
approach 


lim  [/?/,]- 

f  —  C 


1 

2  k*t2 


e 


-2y< 


(115) 


lim[d>]  =  v/2y/k5cos(k(j:-x0)).  (116) 

However,  by  this  time  the  growth  would  have  proceeded  so  far 
that  the  asymptotic  expansions  (1 5)— (21 )  and  the  approximation 
scheme  used  to  derive  (114)  would  no  longer  be  valid.  Neverthe¬ 
less,  the  onset  of  saturation  is  predicted  by  these  considerations, 
and  therefore  it  is  not  valid  to  assume  that  the  linear  Jeans 
instability  is  the  dominant  one,  even  if  its  initial  growth  rate  is 
larger  than  that  of  the  nonlinear  instability.  The  modes  with  k  -»0 
appear  to  reach  the  largest  amplitudes  before  saturating. 


S.  Conclusions 

According  to  Eqs.  (84H86)  and  (102H104),  self-focusing  occurs 
when  w2—0  along  the  acoustic  section.  “Longitudinal"  modula- 
tional  instability  can  occur  on  the  principal  sections  for  small  to 
moderate  values  of  w3,  whereas  the  crossover  regions  (for  the 
case  (n0/,V0)<l)  are  modulationally  stable.  “Transverse”  modu- 
lational  instability  can  occur  at  sufficiently  high  frequencies 
along  both  sections.  For  certain  values  of  the  parameters,  it  is 
also  possible  for  self-focusing  to  occur  along  the  acoustic  section 
at  a  wavenumber  for  which  both  the  acoustic  and  principal 
modes  are  stable.  This  means  that  in  an  expanding  universe, 
where  the  Jeans  length  is  in  fact  slowly  decreasing,  a  wave  of  fixed 
k  will  be  subject  to  modulational  instability  and  self-focusing 
before  the  linear  Jeans  instability  sets  in. 

It  is  of  interest  to  estimate  the  number  of  soliions  of  self- 
focusing  singularities  that  might  be  expected  to  occur  for  a  given 
set  of  initial  conditions.  The  evolution  of  a  soliton  from  some 
initial  data  typically  requires  the  initial  pulse  envelope  to  have 
some  region  over  which  the  phase  of  the  pulse  is  fairly  constant 
and  which  contains  a  normalized  area  of  at  least  n/2.  Over  such  a 
region,  the  typical  number  of  solitons  N ,  is  the  normalized  area 
divided  by  x.  Our  estimate  thus  becomes 

N,*{W/D)'13^  0/io].  017) 

or 

where  /j0  is  roughly  the  size  of  the  envelope  measured  in  the 
original  coordinates,  which  is  related  to  the  degree  of  coherence 
of  the  initial  conditions  and  depends  on  some  correlation  length. 
Along  both  branches  for  large  w,  N,  diverges  like  u>*.  Along  the 
acoustic  section  for  w  —  0.  Nt  diverges  like  to"3,  along  the 
principal  section  like  to" \  For  moderate  values  of  to,  however, 
soliton  formation  is  unlikely  unless  /i0  g  Oft"’]. 


<ArW~ 


toJ 


Wk*c2 


(118) 


whereas  the  typical  time  for  the  dispersion  to  have  an  effect  is 
roughly 


(Arfe 


4*  Vo 

Dw 


(119) 


Along  either  branch  we  find,  for  to  -»  oc, 

(ArW~£'V-7  (120) 

(Arjc-^oW-  (121) 

Along  the  acoustic  section  for  to3  -*  0 

(At)nL~w£'2  (122) 

(At)o~/io“‘3.  (123) 

while  along  the  principal  section  for  small  positive  frequencies 
(Af)*L~CD£'2  (124) 

(Affc-MoW3-  ('25) 


Note  that  it  is  possible  to  achieve  short  nonlinear  growth  time 
scales  for  both  sufficiently  small  and  sufficiently  large  frequencies. 
Since  the  number  of  solitons  estimated  above  also  grows  at  large 
and  small  frequencies,  it  seems  likely  that  modulational  instabi¬ 
lity  and  self-focusing  will  play  a  role  in  a  more  realistic  theory 
which  includes  the  background  time-dependence.  For  moderate 
frequencies,  growth  times  are  long  compared  to  the  time  scale  for 
the  expansion  of  the  universe.  However,  even  in  this  case  it  is 
possible  that  (as  in  the  case  of  the  linear  Jeans  instability)  the 
qualitative  predictions  of  our  “frozen”  background  will  shed 
some  light  on  what  is  to  be  expected  in  a  realistic  background. 
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Appendix:  the  one-fluid  limit 

In  Kates  (1986),  the  two-dimensional  NLS  equation  was  first 
derived  for  a  one-fluid  model.  Unfortunately,  the  coefficients 
given  there  were  incorrect.  Eqs.  (34H38)  of  that  paper  should 


read 

C  =  2  ~^c3k2 

(34) 

D  s  ^[4c,3fc3-3] 

3  k 

(35) 

1  1  , 

£,-2P  +  5f' 

(36) 

-  2icoa, - V  au  +  c2art  +  Qa2a •  *  0 

or 

(37) 
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Q*-o>*--<o>--  =  -c:k*-lc?k2  +  4  (38) 

(Eqs.  (35)  and  (38)  were  correct  as  given  in  the  original  paper) 
Evidently,  the  coefficients  never  take  the  same  sign,  and  self- 
focusing  therefore  does  not  occur.  However,  modulational  inst¬ 
ability  occurs  for  any  m.  Q  changes  sign  at  around  kNL  =  1.336  k,, 
as  reported  in  the  original  paper.  Below  this  critical 
wavenumber,  “longitudinal"  solitons  are  permitted.  Above  kNL, 
“transverse"  solitons  are  permitted. 
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1.  Nonlinear  integrable  lattices  (i.e.  integrable  partial  difference  equations,cf  e.g.  [1,  3]) 
are  of  fundamental  importance  for  the  study  of  classically  integrable  systems.  They  are 
generic  in  the  sense  that  their  various  continuous  limits  give  rise  to  the  hierarchies  of  inte¬ 
grable  PDE’s  [4].  Furthermore,  their  study  opens  up  some  new  points  of  view  on  classical 
integrability  in  general  [5].  In  this  note  we  report  on  another  application  of  such  systems. 
In  fact,  we  will  show  how  these  lattices  give  rise  to  nonlinear  integrable  mappings.  Such 
mappings  are  of  interest  for  the  investigation  of  various  aspects  of  dynamical  systems  (e.g. 
to  study  bifurcations,  transition  to  chaos,  perturbation  techniques,  cf.  e.g.  [7]). 

Probably  the  oldest  nonlinear  integrable  mapping  is  the  elliptic  billiard  due  to  C.G.J.  Ja¬ 
cobi  [6].  More  recently  E.M.  McMillan  found  a  four-parameter  family  of  rational  mappings 
of  the  plane,  together  with  their  invariants  [8].  An  eighteen-parameter  family  generalizing 
the  one  of  McMillan  was  presented  by  G.R.W  Quispel  et  al.  in  [9].  Moreover,  a  connection 
with  soliton  equations  of  differential-difference  type  was  established,  cf.  also  [10].  However, 
a  spectral  interpretation  of  the  integrability  of  these  mappings  on  the  basis  of  a  Lax  pair 
was  lacking. 

In  this  note  we  take  a  different  point  of  view  from  the  one  expounded  in  [9,  10]  by 
considering  integrable  lattices  rather  than  differential-  difference  equations  as  a  starting 
point.  This  is  convenient,  because  it  allows  us  to  obtain  the  mappings  in  a  more  natural 
way  than  before,  namely  not  as  special  reductions,  but  from  the  consideration  of  an  initial 
value  problem  on  a  two-dimensional  planar  lattice.  Furthermore,  we  shall  show  that  these 
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mappings  do  indeed  carry  a  spectral  interpretation,  and  that  the  invariants  can  be  system¬ 
atically  constructed  from  a  Lax  pair. 

2.  We  shall  use  as  a  prime  example  for  the  exposition  of  our  ideas  the  lattice  KdV  equation. 
This  equation  (  as  well  as  lattice  analogs  of  other  integrable  PDE’s  )  was  obtained  in  [1] 
using  a  discrete  version  of  the  direct  linearization  method  introduced  in  [2]. 

The  equation  reads 


(p-g  +  u-u)(p+g-ii+u)  =  p7-q2  .  (1) 

In  (  1)  u  =  u(n,m)  is  the  dynamical  variable  at  the  lattice  site  (n,  m),  n,  m  €  Z,  the  *  and" 
are  shorthand  notations  for  translation  on  the  lattice,  i.e.  ti  =  ti(n  +  1,  m),  ti  =  u(n,  m  +  1) 
and  p  and  q  are  the  lattice  parameters  p,  q  e  C.  Eq.  (  1)  arises  as  the  compatibility 
condition  of  a  pair  of  linear  problems  (Lax  pair)  defining  the  shifts  (translations)  of  an 
eigenfunction  (  k  being  the  spectral  parameter  )  in  the  n-  and  m-directions, 

(p  -  *)* fc  =  Lk-*k,(q-  k)*k  =  Mk  ■  *k  ,  (2) 


where  Lk  is  given  by 

p  -  u  1 

k7  -  p2  +  (p  +  tz)(p-  u)  p+u 

and  where  Mk  is  given  by  a  similar  matrix  obtained  from  (  3)  by  making  the  replacements 
p  —*  q  and  — ►  *. 

Let  us  now  consider  an  initial  value  problem  for  (  1)  on  the  lattice.  One  way  of  doing 
this  is  to  assign  initial  data  on  a  "staircase”  as  in  Fig.  1.  From  the  fact  that  eq.  (  1) 
involves  only  the  four  variables  situated  on  the  four  lattice  sites  around  a  simple  plaquette, 
it  follows  that  the  information  on  these  staircases  evolves  diagonally  through  the  lattice 
along  ’’parallel”  staircases.  Hence,  the  initial  value  problem  is  well-defined. 

Consider, now, the  case  of  periodic  initial  values  along  the  staircases.  The  simplest  non¬ 
trivial  example  is  drawn  in  Fig.  1.  where  we  have  period  2  initial  data  on  two  diagonals, 
a,b,  c,  d  denoting  the  different  initial  values  for  u.  By  applying  the  lattice  equation  we  can 
calculate  the  data  on  all  the  diagonals  where  the  (multiple)  primes  denote  the  various  itera¬ 
tions  of  this  procedure.  One  way  of  doing  this  is  by  regarding  the  first  iteration  as  a  vertical 
shift  on  the  lattice,  thereby  ’’updating”  the  values  of  a, 6, c, das  follows:  b'  =  c,  d'  =  a,  and 
a'  and  c'  are  calculated  in  terms  of  a,6,c,  d  using  (1). 
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Fig  1.  Periodic  configuration  of  initial  data  on  the  lattice  and  its  iteration. 

However, a  reduction  to  a  2-dimensional  one  is  readily  obtained  by  considering  the  dif¬ 
ferences  along  the  diagonalsiAT  =  a  —  c,Y  =  b  -  d  .The  mapping  reads 

y  =  -x ,  X'  =  Y  +  -cr~2  >  (4) 

where  6  =  p-q  and  e  =  p+g.  Eq.  (  4)  is  the  original  McMillan  map  [8]  which  was  recovered 
as  a  reduction  of  the  discrete  Nonlinear  Schrodinger  equation  in  [10].  This  map  arises  as 
the  compatibility  condition  of  a  Lax  pair  that  can  be  constructed  from  (  2)  as  follows  . 

We  introduce  an  ordered  product  of  the  Lax  matrices  (  3)  along  the  staircase.  Because 
of  the  periodicity  it  is  sufficient  to  consider  a  product  of  only  four  of  them,  i.e. 

C.  ~  hdad  •  Ldc  *  '  ^ba  (^) 

where  Lba  and  Ldc  denote  the  matrices  L  of  (  3)  in  which  we  substitute  u  —>  a,u  -*  b 
respectively  u  — >  c,u  —>  d,  and  similarly  in  M  we  substitute  u  — ►  6,  ti  — ►  c  and  u  — ♦  d,ii  — ►  a 
to  obtain  Mcb  and  Mad  respectively.  The  matrix  C  can  be  regarded  as  a  monodromy  matrix 
with  invariant  spectrum  since 

£  =Mala-C-M~l  .  (6) 

In  particular  we  have  that  tr(C)  is  invariant  under  the  mapping,  and  this  indeed  gives  us 
the  McMillan  invariant  [8] 

X  =  (c2  -  X 2)  (e2  -  y2)  +  2 (6XY  ,  (7) 
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which  can  be  parametrized  in  terms  of  Jacobi  elliptic  functions. 


3.  We  have  shown  how  to  obtain  in  a  natural  way  an  integrable  mapping  from  the  lat¬ 
tice  KdV  equation  (1).  Of  course  this  is  only  a  simple  example,  and  there  are  various 
ways  to  generalize  our  procedure.  First  of  all,  we  can  obtain  higher-dimensional  mappings 
related  to  the  lattice  KdV  by  considering  higher  periods  in  the  initial  data.  In  the  case 
of  period  3  for  example  this  yields  a  4-dimensional  mapping  with  two  invariants  that  can 
be  obtained  using  a  monodromy  matrix  as  above.  Another  generalization  is  to  consider 
other  types  of  staircases,  generally  some  discrete  curves  on  the  lattice.  Apart  from  the 
lattice  KdV  one  might  consider  other  existing  lattice  equations  as  starting  point,  such  as 
the  lattice  MKdV,  the  discrete-time  Toda  equation,  the  lattice  BSQ  and  MBSQ  equations 
[11]  and  so  on.  This  work  is  in  progress  [12].  It  would  be  of  interest  to  investigate  the 
canonical  structure  of  the  mappings  and  of  the  lattices  from  which  they  are  obtained.  For 
related  work  towards  this  direction  see  e.g.  [13]. 
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